!**********************************************************************************************************************************
! This file is part of FEAMooring.
!
! Copyright (C) 2014-2016 NREL/TAMU
!
! Licensed under the Apache License, Version 2.0 (the "License");
! you may not use this file except in compliance with the License.
! You may obtain a copy of the License at
!
!     http://www.apache.org/licenses/LICENSE-2.0
!
! Unless required by applicable law or agreed to in writing, software
! distributed under the License is distributed on an "AS IS" BASIS,
! WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
! See the License for the specific language governing permissions and
! limitations under the License.
!
!**********************************************************************************************************************************
    MODULE FEAMooring

    USE FEAMooring_Types
    USE NWTC_Library

    IMPLICIT NONE

    PRIVATE

   TYPE(ProgDesc), PARAMETER            :: FEAM_Ver = ProgDesc( 'FEAMooring', '', '' )

! ==================================================================================================="
! NOTE: The following lines of code were generated by a Matlab script called "Write_ChckOutLst.m"
!      using the parameters listed in the "OutListParameters.xlsx" Excel file. Any changes to these 
!      lines should be modified in the Matlab script and/or Excel worksheet as necessary. 
! ==================================================================================================="
! This code was generated by Write_ChckOutLst.m at 09-Dec-2014 14:03:37.


     ! Indices for computing output channels:
     ! NOTES: 
     !    (1) These parameters are in the order stored in "OutListParameters.xlsx"
     !    (2) Array y%AllOuts() must be dimensioned to the value of the largest output parameter

     !  Time: 

   INTEGER(IntKi), PARAMETER      :: Time      =  0


     ! Category:

   INTEGER(IntKi), PARAMETER      :: FairT1    =  1
   INTEGER(IntKi), PARAMETER      :: FairT2    =  2
   INTEGER(IntKi), PARAMETER      :: FairT3    =  3
   INTEGER(IntKi), PARAMETER      :: FairT4    =  4
   INTEGER(IntKi), PARAMETER      :: FairT5    =  5
   INTEGER(IntKi), PARAMETER      :: FairT6    =  6
   INTEGER(IntKi), PARAMETER      :: FairT7    =  7
   INTEGER(IntKi), PARAMETER      :: FairT8    =  8
   INTEGER(IntKi), PARAMETER      :: FairT9    =  9
   INTEGER(IntKi), PARAMETER      :: FairT10   = 10
   INTEGER(IntKi), PARAMETER      :: AnchT1    = 11
   INTEGER(IntKi), PARAMETER      :: AnchT2    = 12
   INTEGER(IntKi), PARAMETER      :: AnchT3    = 13
   INTEGER(IntKi), PARAMETER      :: AnchT4    = 14
   INTEGER(IntKi), PARAMETER      :: AnchT5    = 15
   INTEGER(IntKi), PARAMETER      :: AnchT6    = 16
   INTEGER(IntKi), PARAMETER      :: AnchT7    = 17
   INTEGER(IntKi), PARAMETER      :: AnchT8    = 18
   INTEGER(IntKi), PARAMETER      :: AnchT9    = 19
   INTEGER(IntKi), PARAMETER      :: AnchT10   = 20

     ! Category:


     ! The maximum number of output channels which can be output by the code.
   INTEGER(IntKi), PARAMETER      :: MaxOutPts = 20

!End of code generated by Matlab script

   INTEGER(IntKi), PARAMETER         :: FairT(10) = (/ FairT1, FairT2, FairT3, FairT4, FairT5, FairT6, FairT7, FairT8, FairT9, FairT10 /)
   INTEGER(IntKi), PARAMETER         :: AnchT(10) = (/ AnchT1, AnchT2, AnchT3, AnchT4, AnchT5, AnchT6, AnchT7, AnchT8, AnchT9, AnchT10 /)


    ! ..... Public Subroutines ...................................................................................................

    PUBLIC :: FEAM_Init                           ! Initialization routine
    PUBLIC :: FEAM_Solve
    PUBLIC :: FEAM_End                            ! Ending routine (includes clean up)
    PUBLIC :: FEAM_UpdateStates                   ! Loose coupling routine for solving for constraint states, integrating
    PUBLIC :: FEAM_CalcOutput                     ! Routine for computing outputs


    CONTAINS
!----------------------------------------------------------------------------------------------------------------------------------
SUBROUTINE FEAM_Init(  InitInp, u, p, x, xd, z, OtherState, y, misc, Interval, InitOut, ErrStat, ErrMsg )
! This routine is called at the start of the simulation to perform initialization steps.
! The parameters are set here and not changed during the simulation.
! The initial states and initial guess for the input are defined.
!..................................................................................................................................

    TYPE(FEAM_InitInputType),       INTENT(IN   )  :: InitInp     ! Input data for initialization routine
    TYPE(FEAM_InputType),           INTENT(  OUT)  :: u           ! An initial guess for the input; input mesh must be defined
    TYPE(FEAM_ParameterType),       INTENT(  OUT)  :: p           ! Parameters
    TYPE(FEAM_ContinuousStateType), INTENT(  OUT)  :: x           ! Initial continuous states
    TYPE(FEAM_DiscreteStateType),   INTENT(  OUT)  :: xd          ! Initial discrete states.
    TYPE(FEAM_ConstraintStateType), INTENT(  OUT)  :: z           ! Initial guess of the constraint states.
    TYPE(FEAM_OtherStateType),      INTENT(  OUT)  :: OtherState  ! Initial other states
    TYPE(FEAM_OutputType),          INTENT(  OUT)  :: y           ! Initial system outputs (outputs are not calculated;
                                                                  !   only the output mesh is initialized)
    TYPE(FEAM_MiscVarType),         INTENT(  OUT)  :: misc        ! Initial misc/optimization variables
    REAL(DbKi),                     INTENT(INOUT)  :: Interval    ! Coupling interval in seconds: the rate that
                                                                  !   (1) FEAM_UpdateStates() is called in loose coupling &
                                                                  !   (2) FEAM_UpdateDiscState() is called in tight coupling.
                                                                  !   Input is the suggested time from the glue code;
                                                                  !   Output is the actual coupling interval that will be used
                                                                  !   by the glue code.
    TYPE(FEAM_InitOutputType),      INTENT(  OUT)  :: InitOut     ! Output for initialization routine
    INTEGER(IntKi),                 INTENT(  OUT)  :: ErrStat     ! Error status of the operation
    CHARACTER(*),                   INTENT(  OUT)  :: ErrMsg      ! Error message if ErrStat /= ErrID_None

    ! local variables

    TYPE(FEAM_InputFile)                           :: InputFileData  ! Data stored in the module's input file
    INTEGER(IntKi)                                 :: ErrStat2       ! temporary Error status of the operation
    CHARACTER(ErrMsgLen)                           :: ErrMsg2        ! temporary Error message if ErrStat /= ErrID_None

    REAL(ReKi)                                     :: Dist           ! A temporary real  value
    REAL(ReKi)                                     :: Curv           ! A temporary real  value
    REAL(ReKi)                                     :: Temp           ! A temporary real  value
    INTEGER(IntKi)                                 :: I              ! Generic loop index
    INTEGER(IntKi)                                 :: J              ! Generic loop index


    ! Initialize variables

    ErrStat = ErrID_None
    ErrMsg  = ""

    ! Initialize the NWTC Subroutine Library

    CALL NWTC_Init( EchoLibVer=.FALSE. )

    ! Display the module information

    CALL DispNVD( FEAM_Ver )


    !............................................................................................      
    ! Read the input file and validate the data
    ! (note p%NumBl and p%RootName must be set first!) 
    !............................................................................................      
    p%RootName = TRIM(InitInp%RootName) ! all of the output file names from this module will end with '.FEAM.*'   

    CALL FEAM_ReadInput( InitInp%InputFile, InputFileData, p%RootName, Interval, InitInp%gravity, InitInp%WtrDens, ErrStat2, ErrMsg2 )
    CALL CheckError( ErrStat2, ErrMsg2 )
    IF (ErrStat >= AbortErrLev) RETURN

    !............................................................................................
    ! Define parameters here:
    !............................................................................................
    CALL FEAM_SetParameters( InputFileData, p, ErrStat2, ErrMsg2 )   
    CALL CheckError( ErrStat2, ErrMsg2 )
    IF (ErrStat >= AbortErrLev) RETURN      

    !............................................................................................
    ! Define initialization-routine output here (we're doing this near the beginning because  
    !   we use LAnch{x|y|z}i and LFair{x|y|z}t  to compute some parameters:
    !............................................................................................

    CALL Init_InitOut(p, InputFileData, InitInp, InitOut, ErrStat2, ErrMsg2)
      CALL CheckError( ErrStat2, ErrMsg2 )
      IF (ErrStat >= AbortErrLev) RETURN
      
      
    !............................................................................................
    ! Define additional parameters based on InitOut data and necessary for state initializion here:
    !............................................................................................
      
    DO I = 1,p%NumLines
        Dist = SQRT( ( InitOut%LFairxt(I) - InitOut%LAnchxi(I) )**2 + &
                     ( InitOut%LFairyt(I) - InitOut%LAnchyi(I) )**2 + &
                     ( InitOut%LFairzt(I) - InitOut%LAnchzi(I) )**2 )

        Curv = 0.0
        DO J = 1,p%NumElems
            Curv = Curv + p%Elength(I)
        ENDDO  

        p%GLUZR(I,1,1) = InitOut%LAnchxi(I)
        p%GLUZR(I,1,3) = InitOut%LAnchyi(I)
        p%GLUZR(I,1,5) = InitOut%LAnchzi(I)

        p%GLUZR(I,1,2) = (InitOut%LFairxt(I) - InitOut%LAnchxi(I))/Dist
        p%GLUZR(I,1,4) = (InitOut%LFairyt(I) - InitOut%LAnchyi(I))/Dist
        p%GLUZR(I,1,6) = (InitOut%LFairzt(I) - InitOut%LAnchzi(I))/Dist

        Temp = 0.0
        DO J = 2,p%NumElems+1
            Temp = Temp+p%Elength(I)
            p%GLUZR(I,J,1) = InitOut%LAnchxi(I)+(InitOut%LFairxt(I)- InitOut%LAnchxi(I))*Temp/Curv
            p%GLUZR(I,J,3) = InitOut%LAnchyi(I)+(InitOut%LFairyt(I)- InitOut%LAnchyi(I))*Temp/Curv
            p%GLUZR(I,J,5) = InitOut%LAnchzi(I)+(InitOut%LFairzt(I)- InitOut%LAnchzi(I))*Temp/Curv

            p%GLUZR(I,J,2) = (InitOut%LFairxt(I)-InitOut%LAnchxi(I))/Dist
            p%GLUZR(I,J,4) = (InitOut%LFairyt(I)-InitOut%LAnchyi(I))/Dist
            p%GLUZR(I,J,6) = (InitOut%LFairzt(I)-InitOut%LAnchzi(I))/Dist
        ENDDO  
    ENDDO  
      
          
    DO I = 1,p%NumLines
        p%GP(I,1) = InitOut%LFairxt(I)
        p%GP(I,2) = InitOut%LFairyt(I)
        p%GP(I,3) = InitOut%LFairzt(I)
    ENDDO  

        ! INPUT BOUNDARY VALUE
    DO I = 1,p%NumLines
        p%Bvp(I,1) = InitOut%LAnchxi(I)
        p%Bvp(I,2) = InitOut%LAnchyi(I)
        p%Bvp(I,3) = InitOut%LAnchzi(I)
    ENDDO  

    !............................................................................................
    ! Define initial system states here:
    !............................................................................................

    CALL Init_MiscVars(p, misc, ErrStat2, ErrMsg2 )    
      CALL CheckError( ErrStat2, ErrMsg2 )
      IF (ErrStat >= AbortErrLev) RETURN
            
    CALL Init_States(p, x, xd, z, OtherState, ErrStat2, ErrMsg2 )    
      CALL CheckError( ErrStat2, ErrMsg2 )
      IF (ErrStat >= AbortErrLev) RETURN
    
    CALL FEAM_Solve(0.0_Dbki, 0, u, p, x, xd, z, OtherState, misc, .false., ErrStat2, ErrMsg2 )
      CALL CheckError( ErrStat2, ErrMsg2 )
      IF (ErrStat >= AbortErrLev) RETURN       
    

    ! Revert back fairlead position for dynamic analysis (and mesh initialization)
    DO I = 1,p%NumLines
        p%GP(I,1) = InitOut%LFairxt(I) - InitInp%PtfmInit(1)
        p%GP(I,2) = InitOut%LFairyt(I) - InitInp%PtfmInit(2)
        p%GP(I,3) = InitOut%LFairzt(I) - InitInp%PtfmInit(3)
    END DO

    
    !............................................................................................
    ! Define initial guess for the system inputs here:
    !............................................................................................
        
    ! Create Line mesh
    CALL CreateLineMesh( u, p, misc, ErrStat, ErrMsg )

    ! Create Platform point mesh
    CALL CreatePtFairleadMesh( u, p, y, ErrStat, ErrMsg ) 
    
       !.........................................................................................
       ! initialization for WaveAcc0, etc, which TAMU had in HydroDyn
       !.........................................................................................
       p%NStepWave = InitInp%NStepWave
       ALLOCATE( p%WaveAcc0 ( 0:p%NStepWave, u%HydroForceLineMesh%NNodes, 3), STAT=ErrStat2 )
         IF (ErrStat2 /= 0) CALL CheckError(ErrID_Fatal, "FEAM_Init:Error allocating p%WaveAcc0z")         
       ALLOCATE( p%WaveTime ( 0:p%NStepWave                                ), STAT=ErrStat2 )
         IF (ErrStat2 /= 0) CALL CheckError(ErrID_Fatal, "FEAM_Init:Error allocating p%WaveTime")
       ALLOCATE( p%WaveVel0 ( 0:p%NStepWave, u%HydroForceLineMesh%NNodes, 3), STAT=ErrStat2 )
         IF (ErrStat2 /= 0) CALL CheckError(ErrID_Fatal, "FEAM_Init:Error allocating p%WaveVel0")
    
       DO i=0,p%NStepWave
          p%WaveTime = i*p%DT  !bjj: just made up some monotonic sequence here
       END DO
    
       p%WaveAcc0 = 0.0_ReKi
       p%WaveVel0 = 0.0_ReKi
        
    
      
    
    !............................................................................................
    ! Define system output initializations (set up mesh) here:
    !............................................................................................
    CALL AllocAry( y%WriteOutput, p%NumOuts, 'WriteOutput', ErrStat2, ErrMsg2 )
    CALL CheckError( ErrStat2, ErrMsg2 )
    IF (ErrStat >= AbortErrLev) RETURN
    y%WriteOutput = 0

    !............................................................................................
    ! If you want to choose your own rate instead of using what the glue code suggests, tell the glue code the rate at which
    !   this module must be called here:
    !............................................................................................

    Interval = p%DT

    ! Clean up the local variable:
    CALL FEAM_DestroyInputFile( InputFileData, ErrStat2, ErrMsg2 )

    RETURN

CONTAINS
    !...............................................................................................................................
    SUBROUTINE CheckError(ErrID,Msg)
    ! This subroutine sets the error message and level and cleans up if the error is >= AbortErrLev
    !...............................................................................................................................

    ! Passed arguments
    INTEGER(IntKi), INTENT(IN) :: ErrID       ! The error identifier (ErrStat)
    CHARACTER(*),   INTENT(IN) :: Msg         ! The error message (ErrMsg)

    INTEGER(IntKi)             :: ErrStat3    ! The error identifier (ErrStat)
    CHARACTER(ErrMsgLen)       :: ErrMsg3     ! The error message (ErrMsg)

    !............................................................................................................................
    ! Set error status/message;
    !............................................................................................................................

    IF ( ErrID /= ErrID_None ) THEN

    IF ( LEN_TRIM(ErrMsg) > 0 ) ErrMsg = TRIM(ErrMsg)//NewLine
    ErrMsg = TRIM(ErrMsg)//'Error in FEAM_Init: '//TRIM(Msg)
    ErrStat = MAX(ErrStat, ErrID)

    !.........................................................................................................................
    ! Clean up if we're going to return on error: close files, deallocate local arrays
    !.........................................................................................................................
    IF ( ErrStat >= AbortErrLev ) THEN
        CALL FEAM_DestroyInputFile(InputFileData, ErrStat3, ErrMsg3 )
    END IF

    END IF

    END SUBROUTINE CheckError 
    !----------------------------------------------------------------------------------------------------------------------------------
END SUBROUTINE FEAM_Init
!----------------------------------------------------------------------------------------------------------------------------------
SUBROUTINE Init_States(p, x, xd, z, OtherState, ErrStat, ErrMsg)

   TYPE(FEAM_ParameterType),       INTENT(IN   )  :: p           ! Parameters
   TYPE(FEAM_ContinuousStateType), INTENT(INOUT)  :: x           ! Initial continuous states
   TYPE(FEAM_DiscreteStateType),   INTENT(INOUT)  :: xd          ! Initial discrete states. 
   TYPE(FEAM_ConstraintStateType), INTENT(INOUT)  :: z           ! Initial guess of the constraint states.
   TYPE(FEAM_OtherStateType),      INTENT(INOUT)  :: OtherState  ! Initial other states

   INTEGER(IntKi),                 INTENT(  OUT)  :: ErrStat     ! Error status of the operation
   CHARACTER(*),                   INTENT(  OUT)  :: ErrMsg      ! Error message if ErrStat /= ErrID_None
   
    ! local variables
   INTEGER(IntKi)                                 :: ErrStat2     ! Error status of the operation
   CHARACTER(ErrMsgLen)                           :: ErrMsg2      ! Error message if ErrStat /= ErrID_None
   CHARACTER(*), PARAMETER                        :: RoutineName = 'Init_States'
   
   INTEGER(IntKi)                                 :: I,  J,  K
   INTEGER(IntKi)                                 :: MM, NN, NE   ! counters
   
   ErrStat = ErrID_None
   ErrMsg = ""
   
      ! continuous states
   CALL AllocAry( x%GLU,  p%NumLines, 8*p%NumElems+7, 'GLU',  ErrStat2, ErrMsg2 ); call SetErrStat(ErrStat2,ErrMsg2,ErrStat,ErrMsg,RoutineName)
   CALL AllocAry( x%GLDU, p%NumLines, 8*p%NumElems+7, 'GLDU', ErrStat2, ErrMsg2 ); call SetErrStat(ErrStat2,ErrMsg2,ErrStat,ErrMsg,RoutineName)
   
   IF (ErrStat >= AbortErrLev) RETURN

         
      ! discrete states
   
    xd%DummyDiscState = 0
   
      ! constraint states
   
   ! these are initialized in FEAM_Solve()
    
      ! other states: 

   CALL AllocAry( OtherState%GLU0,    p%NumLines,8*p%NumElems+7,            'GLU0',    ErrStat2,ErrMsg2); call SetErrStat(ErrStat2,ErrMsg2,ErrStat,ErrMsg,RoutineName)
   CALL AllocAry( OtherState%GLDDU,   p%NumLines,8*p%NumElems+7,            'GLDDU',   ErrStat2,ErrMsg2); call SetErrStat(ErrStat2,ErrMsg2,ErrStat,ErrMsg,RoutineName)
   CALL AllocAry( OtherState%GFORC0,  p%NumLines,p%NumElems,p%NBAND,        'GFORC0',  ErrStat2,ErrMsg2); call SetErrStat(ErrStat2,ErrMsg2,ErrStat,ErrMsg,RoutineName)
   CALL AllocAry( OtherState%GMASS0,  p%NumLines,p%NumElems,p%NBAND,p%NBAND,'GMASS0',  ErrStat2,ErrMsg2); call SetErrStat(ErrStat2,ErrMsg2,ErrStat,ErrMsg,RoutineName)
   CALL AllocAry( OtherState%FAST_FPA,p%NumLines,p%NDIM,                    'FAST_FPA',ErrStat2,ErrMsg2); call SetErrStat(ErrStat2,ErrMsg2,ErrStat,ErrMsg,RoutineName)
   CALL AllocAry( OtherState%FAST_RP, p%NumLines,p%NDIM,                    'FAST_RP', ErrStat2,ErrMsg2); call SetErrStat(ErrStat2,ErrMsg2,ErrStat,ErrMsg,RoutineName)   
   
   if (ErrStat >= AbortErrLev) return

   OtherState%BottomTouch = .FALSE.  
    
   OtherState%GFORC0      = 0.0_ReKi
   OtherState%GMASS0      = 0.0_ReKi
    
   OtherState%FAST_FPA    = p%GP    
   OtherState%FAST_RP     = 0.0 
    
   OtherState%INCR        = 0
    
   OtherState%RSDF        = 0.0_ReKi
    
   OtherState%FORC0       = 0.0_ReKi
   OtherState%EMAS0       = 0.0_ReKi

       
   OtherState%GLU0   = 0.0_ReKi
   OtherState%GLDDU  = 0.0_ReKi
   x%GLU  = 0.0
   x%GLDU = 0.0   

   DO I = 1,p%NumLines

      !SET THE INITIAL VALUE TO GLOBE VECTOR GLU
      DO NN = 1,p%NumElems+1
         DO J = 1,p%NDIM
            DO K = 1,2
               MM = p%IDOF(J,K)+(NN-1)*p%NHBD
                        x%GLU(  I,MM) = p%GLUZR(I,NN,p%IDOF(J,K))
               OtherState%GLU0( I,MM) = x%GLU(I,MM)
                        x%GLDU( I,MM) = 0.
               OtherState%GLDDU(I,MM) = 0.
            END DO
         END DO
      END DO

      DO NE = 1,p%NumElems
         DO J = 1,3
               MM = p%JDOF(J)+(NE-1)*p%NHBD
               x%GLU(I,MM)  = p%GTZER(I,J+2*(NE-1))
               OtherState%GLU0(I,MM) = x%GLU(I,MM)
         END DO
      END DO

   ENDDO      
   
END SUBROUTINE Init_States
!----------------------------------------------------------------------------------------------------------------------------------
SUBROUTINE Init_miscVars(p,misc,ErrStat,ErrMsg)
   TYPE(FEAM_ParameterType),       INTENT(IN   )  :: p           ! Parameters
   TYPE(FEAM_MiscVarType),         INTENT(INOUT)  :: misc        ! misc/optimization variables

   INTEGER(IntKi),                 INTENT(  OUT)  :: ErrStat     ! Error status of the operation
   CHARACTER(*),                   INTENT(  OUT)  :: ErrMsg      ! Error message if ErrStat /= ErrID_None
   
    ! local variables
   INTEGER(IntKi)                                 :: ErrStat2     ! Error status of the operation
   CHARACTER(ErrMsgLen)                           :: ErrMsg2      ! Error message if ErrStat /= ErrID_None
   CHARACTER(*), PARAMETER                        :: RoutineName = 'Init_miscVars'
   
   ErrStat = ErrID_None
   ErrMsg = ""
   
   CALL AllocAry(misc%Line_Coordinate,p%NumLines, p%NumElems+1, p%NDIM,'Line_Coordinate',ErrStat2,ErrMsg2); call SetErrStat(ErrStat2,ErrMsg2,ErrStat,ErrMsg,RoutineName)
   CALL AllocAry(misc%Line_Tangent,   p%NumLines, p%NumElems+1, p%NDIM,'Line_Tangent',   ErrStat2,ErrMsg2); call SetErrStat(ErrStat2,ErrMsg2,ErrStat,ErrMsg,RoutineName)
   CALL AllocAry(misc%GLF,            p%NumLines, 8*p%NumElems+7,        'GLF',          ErrStat2,ErrMsg2); call SetErrStat(ErrStat2,ErrMsg2,ErrStat,ErrMsg,RoutineName)
   CALL AllocAry(misc%GLK,            p%NumLines, 8*p%NumElems+7,p%NBAND,'GLK',          ErrStat2,ErrMsg2); call SetErrStat(ErrStat2,ErrMsg2,ErrStat,ErrMsg,RoutineName)
   CALL AllocAry(misc%FAST_FP,        p%NumLines, p%NDIM,                'FAST_FP',      ErrStat2,ErrMsg2); call SetErrStat(ErrStat2,ErrMsg2,ErrStat,ErrMsg,RoutineName) ! holds absolute position of input mesh
   CALL AllocAry(misc%FAIR_T,         p%NumLines,                        'FAIR_T',       ErrStat2,ErrMsg2); call SetErrStat(ErrStat2,ErrMsg2,ErrStat,ErrMsg,RoutineName)
   CALL AllocAry(misc%ANCH_T,         p%NumLines,                        'ANCH_T',       ErrStat2,ErrMsg2); call SetErrStat(ErrStat2,ErrMsg2,ErrStat,ErrMsg,RoutineName)
   CALL AllocAry(misc%FAIR_ANG,       p%NumLines, p%NDIM,                'FAIR_ANG',     ErrStat2,ErrMsg2); call SetErrStat(ErrStat2,ErrMsg2,ErrStat,ErrMsg,RoutineName)
   CALL AllocAry(misc%ANCH_ANG,       p%NumLines, p%NDIM,                'ANCH_ANG',     ErrStat2,ErrMsg2); call SetErrStat(ErrStat2,ErrMsg2,ErrStat,ErrMsg,RoutineName)
   CALL AllocAry(misc%F_Lines,        p%NumLines, p%NDIM,                'F_Lines',      ErrStat2,ErrMsg2); call SetErrStat(ErrStat2,ErrMsg2,ErrStat,ErrMsg,RoutineName)
   
   if (ErrStat >= AbortErrLev) return
      
   misc%LastIndWave     = 1
   
   
END SUBROUTINE Init_miscVars
!----------------------------------------------------------------------------------------------------------------------------------
SUBROUTINE Init_InitOut(p,InputFileData,InitInp, InitOut,ErrStat,ErrMsg)
   TYPE(FEAM_ParameterType),       INTENT(IN   )  :: p               ! Parameters
   TYPE(FEAM_InputFile),           INTENT(IN   )  :: InputFileData   ! Data stored in the module's input file
   TYPE(FEAM_InitInputType),       INTENT(IN   )  :: InitInp         ! Input data for initialization routine
   TYPE(FEAM_InitOutputType),      INTENT(INOUT)  :: InitOut         ! Output for initialization routine   
   
   INTEGER(IntKi),                 INTENT(  OUT)  :: ErrStat         ! Error status of the operation
   CHARACTER(*),                   INTENT(  OUT)  :: ErrMsg          ! Error message if ErrStat /= ErrID_None
   
    ! local variables
   INTEGER(IntKi)                                 :: ErrStat2        ! Error status of the operation
   CHARACTER(ErrMsgLen)                           :: ErrMsg2         ! Error message if ErrStat /= ErrID_None
   CHARACTER(*), PARAMETER                        :: RoutineName = 'Init_InitOut'
   
   integer(intKi)                                 :: I               ! loop counter
   
   ErrStat = ErrID_None
   ErrMsg = ""
   
   
   InitOut%Ver = FEAM_Ver     
   
   CALL AllocAry(InitOut%WriteOutputHdr,p%NumOuts,'WriteOutputHdr',ErrStat2,ErrMsg2); call SetErrStat(ErrStat2,ErrMsg2,ErrStat,ErrMsg,RoutineName)
   CALL AllocAry(InitOut%WriteOutputUnt,p%NumOuts,'WriteOutputUnt',ErrStat2,ErrMsg2); call SetErrStat(ErrStat2,ErrMsg2,ErrStat,ErrMsg,RoutineName)

   CALL AllocAry(InitOut%LAnchxi,  p%NumLines, 'LAnchxi',  ErrStat2, ErrMsg2); call SetErrStat(ErrStat2,ErrMsg2,ErrStat,ErrMsg,RoutineName)
   CALL AllocAry(InitOut%LAnchyi,  p%NumLines, 'LAnchyi',  ErrStat2, ErrMsg2); call SetErrStat(ErrStat2,ErrMsg2,ErrStat,ErrMsg,RoutineName)
   CALL AllocAry(InitOut%LAnchzi,  p%NumLines, 'LAnchzi',  ErrStat2, ErrMsg2); call SetErrStat(ErrStat2,ErrMsg2,ErrStat,ErrMsg,RoutineName)
   CALL AllocAry(InitOut%LFairxt,  p%NumLines, 'LFairxt',  ErrStat2, ErrMsg2); call SetErrStat(ErrStat2,ErrMsg2,ErrStat,ErrMsg,RoutineName)
   CALL AllocAry(InitOut%LFairyt,  p%NumLines, 'LFairyt',  ErrStat2, ErrMsg2); call SetErrStat(ErrStat2,ErrMsg2,ErrStat,ErrMsg,RoutineName)
   CALL AllocAry(InitOut%LFairzt,  p%NumLines, 'LFairzt',  ErrStat2, ErrMsg2); call SetErrStat(ErrStat2,ErrMsg2,ErrStat,ErrMsg,RoutineName)
   
   if (ErrStat >= AbortErrLev) return
   
   
   DO I = 1,p%NumLines
      InitOut%LAnchxi(I) =  InputFileData%LRadAnch(I) *COS(InputFileData%LAngAnch(I))
      InitOut%LAnchyi(I) =  InputFileData%LRadAnch(I) *SIN(InputFileData%LAngAnch(I))                       ! Convert the radius, azimuth angle, and depth of the current anchor   to x-, y-, and z-coordinates into the inertial frame coordinate system
      InitOut%LAnchzi(I) = -InputFileData%LDpthAnch(I)                                                      !
      InitOut%LFairxt(I) =  InputFileData%LRadFair(I) *COS(InputFileData%LAngFair(I)) + InitInp%PtfmInit(1) ! Add platform initial position to obtain static fairlead offset at initial time
      InitOut%LFairyt(I) =  InputFileData%LRadFair(I) *SIN(InputFileData%LAngFair(I)) + InitInp%PtfmInit(2) ! Add platform initial position to obtain static fairlead offset at initial time
      InitOut%LFairzt(I) = -InputFileData%LDrftFair(I)                                + InitInp%PtfmInit(3) ! Add platform initial position to obtain static fairlead offset at initial time
   END DO
   
   do i=1,p%NumOuts    
      InitOut%WriteOutputHdr(i) = p%OutParam(i)%Name
      InitOut%WriteOutputUnt(i) = p%OutParam(i)%Units     
   end do
   
   
END SUBROUTINE Init_InitOut
!----------------------------------------------------------------------------------------------------------------------------------
SUBROUTINE FEAM_Solve( t, n, u, p, x, xd, z, OtherState, misc, Dyn, ErrStat, ErrMsg )

   TYPE(FEAM_InputType),           INTENT(INOUT)  :: u           ! An initial guess for the input; input mesh must be defined
   TYPE(FEAM_ParameterType),       INTENT(IN   )  :: p           ! Parameters
   TYPE(FEAM_ContinuousStateType), INTENT(INOUT)  :: x           ! Initial continuous states
   TYPE(FEAM_DiscreteStateType),   INTENT(INOUT)  :: xd          ! Initial discrete states. 
   TYPE(FEAM_ConstraintStateType), INTENT(INOUT)  :: z           ! Initial guess of the constraint states.
   TYPE(FEAM_OtherStateType),      INTENT(INOUT)  :: OtherState  ! other states
   TYPE(FEAM_MiscVarType),         INTENT(INOUT)  :: misc        ! misc/optimization variables

   LOGICAL       ,                 INTENT(IN   )  :: Dyn         ! flag for dynamic or static analysis
   INTEGER(IntKi),                 INTENT(  OUT)  :: ErrStat     ! Error status of the operation
   CHARACTER(*),                   INTENT(  OUT)  :: ErrMsg      ! Error message if ErrStat /= ErrID_None

   ! Passed variables:
!   INTEGER(IntKi)                                 :: ErrStat2       ! temporary Error status of the operation
!   CHARACTER(ErrMsgLen)                           :: ErrMsg2        ! temporary Error message if ErrStat /= ErrID_None
   REAL(DbKi),                     INTENT(IN)     :: t              ! Current time
   INTEGER(IntKi),                 INTENT(IN)     :: n              ! Current step
   REAL(ReKi)                   :: RHSR0(6)                         ! Temporary variable (old right hand side R)
   REAL(ReKi)                   :: C0,C00                           ! Temporary variable
   REAL(ReKi)                   :: VP                               ! Temporary variable
   REAL(ReKi)                   :: SUM1,SUM2,SUM3,SUMA,SUMB,SUMC    ! Temporary variable
   REAL(ReKi)                   :: ERR1,ERR2,ERR3                   ! Temporary variable
   REAL(ReKi)                   :: IVERT,TEMP                       ! Temporary variable
   
    
   INTEGER(IntKi)                  :: MM,NN,LEG,NE                  ! Generic index
   INTEGER(IntKi)                  :: I,J,K,L                       ! Generic index
!    INTEGER(IntKi)                  :: Sttus                         ! Status returned by an attempted allocation or READ.
    
   LOGICAL  :: Iteration1      ! Static iteration flag [-]
   LOGICAL  :: Iteration2      ! Static iteration flag [-]

    
    ! Initialize variables

   ErrStat = ErrID_None
   ErrMsg  = ""

    !  SET A COEFFICIENT TO DOUBLE BOUNDARY FORCE IN TIME DOMAIN
    !  AND DN IN STATISTICS CALCULATION

   IF(DYN) THEN
      C0  = 2.0
      C00 = 2.0/p%DT
   ELSE
      C0  = 1.0
      C00 = 0.0
   ENDIF

   Iteration1 = .FALSE.  
   Iteration2 = .FALSE.  
        
   DO WHILE ((.NOT. Iteration1).OR.(.NOT. Iteration2))

       DO LEG = 1, p%NumLines ! Start LEG Loop     

          !  INITAILIZE FORCE VECTOR AND STIFFNESS MATRIX 
          DO I = 1,p%NEQ(LEG)
              misc%GLF(LEG,I) = 0.0
              DO J = 1,p%NBAND
                  misc%GLK(LEG,I,J)=0.0
              ENDDO
          ENDDO

          DO NE = 1, p%NumElems ! Start Element Loop

             ! 1. INPUT DISPLACEMENT
             DO I = 1,p%NDIM
                 DO J = 1,4
                     MM=p%IDOF(I,J)+(NE-1)*p%NHBD
                     misc%U ( I,J)=          x%GLU( LEG,MM)
                     misc%U0( I,J)=OtherState%GLU0( LEG,MM)
                     misc%DU( I,J)=         x%GLDU( LEG,MM)
                     misc%DDU(I,J)=OtherState%GLDDU(LEG,MM)
                 ENDDO
             ENDDO

             ! 2. INPUT LAGRANGIAN MULTIPLIER
             DO I = 1,3
                 MM=p%JDOF(I)+(NE-1)*p%NHBD
                 z%TSN(I)=x%GLU(LEG,MM)
                 z%TZER(I)=p%GTZER(LEG,I+2*(NE-1))
             ENDDO

             ! 3. INPUT EXTERNAL FORCE AND MASS OF PREVIOUS TIME STEP
             IF(OtherState%INCR.NE.0) THEN
                 DO I = 1,p%NDIM
                     DO J = 1,4
                         MM=p%IDOF(I,J)
                         OtherState%FORC0(MM)=OtherState%GFORC0(LEG,NE,MM)
                         DO K = 1,p%NDIM
                             DO L = 1,4
                                 NN=p%IDOF(K,L)
                                 OtherState%EMAS0(MM,NN)=OtherState%GMASS0(LEG,NE,MM,NN)
                             ENDDO
                         ENDDO
                     ENDDO
                 ENDDO
             !ELSE ! IF OtherState%INCR == 0
             !    BJJ 7-Jan-2016: Unless I'm missing something, I don't see how OtherState%FORC0 and OtherState%EMAS0 get 
             !    set in this case... and both values are used in the call to Element(). 
             ENDIF 
             
       ! 4. AND OTHER INFORMATION NEEDED IN ELEMNT
    
             CALL Element(t, n, LEG, NE, u, p, x, xd, z, OtherState, misc, ErrStat, ErrMsg)

                !  SAVE THE EXTERNAL FORCE AND MASS FOR NEXT TIME STEP
                !  AND GET THE NODAL INTERACTION FORCE (STATIC & DYNAMIC) 

             DO I = 1,p%NDIM
                 DO J = 1,4
                     MM=p%IDOF(I,J)
                     OtherState%GFORC0(LEG,NE,MM)=misc%FORCE(MM)
                     DO K = 1,p%NDIM
                         DO L = 1,4
                             NN=p%IDOF(K,L)
                             OtherState%GMASS0(LEG,NE,MM,NN)=misc%EMASS(MM,NN)
                         ENDDO
                     ENDDO
                 ENDDO
             ENDDO  

             IF(NE.EQ.p%NumElems)THEN
                 DO I = 1,p%NDIM
                     misc%SLIN (I) = p%GSL(LEG,2,I)
                     misc%FP   (I) = p%GP(LEG,I)
                     misc%R    (I) = misc%U(I,3)
                     misc%RP   (I) = misc%U(I,4)/p%Elength(LEG)
                 ENDDO

                 DO I=1,p%NDIM
                     OtherState%FAST_RP(LEG,I)  = misc%RP(I)
                 ENDDO  

                 ! First run for BXA             
                 IF(DYN) THEN
                     DO I = 1,p%NDIM
                         misc%FP(I) = OtherState%FAST_FPA(LEG,I)  ! bjj: OtherState%FAST_FPA are values at t (couldn't we just replace them with call to extrap-interp?)
                     ENDDO
                 ENDIF                

                 CALL Couple (p, x, OtherState, misc, ErrStat, ErrMsg)

                 DO I = 1,2*p%NDIM
                     RHSR0(I)=misc%RHSR(I)
                 END DO

                 ! Second run for BX
                 IF(DYN) THEN
                     DO I = 1,p%NDIM
                         misc%FP(I) = misc%FAST_FP(LEG,I)    ! bjj: misc%FAST_FPA are values at t+dt 
                     ENDDO
                 ENDIF 

                 CALL Couple (p, x, OtherState, misc, ErrStat, ErrMsg)

                 DO I = 1,2*p%NDIM
                     MM=p%NHBD+I
                     OtherState%RSDF(MM)=OtherState%RSDF(MM)+misc%RHSR(I)
                     IF(DYN) OtherState%RSDF(MM)=OtherState%RSDF(MM)+RHSR0(I)
                     DO J = 1,2*p%NDIM
                         NN=p%NHBD+J
                         misc%ESTIF(MM,NN)=misc%ESTIF(MM,NN)+misc%STIFR(I,J)
                     ENDDO
                 ENDDO           
             ENDIF

             !  ASSEMBLE THE STIFF MATRIX AND FORCE VECTOR   
             DO I = 1,p%NBAND
                 MM=I+(NE-1)*p%NHBD
                 misc%GLF(LEG,MM)=misc%GLF(LEG,MM)+OtherState%RSDF(I)
                 DO J = I,p%NBAND
                     NN=J-I+1
                     misc%GLK(LEG,MM,NN)=misc%GLK(LEG,MM,NN)+misc%ESTIF(I,J)
                 ENDDO
             ENDDO
         ENDDO ! End Element loop

          !  APPLY ESSENTIAL BOUNDARY CONDITIONS TO GLOBAL MATRIX
          DO L = 1,p%NDIM
              MM=p%IDOF(L,1)
              VP=p%Bvp(LEG,L)-x%GLU(LEG,MM)         
                 IF(OtherState%INCR.NE.0.OR. .NOT. DYN) VP=0.
              I=MM-p%NBAND
              DO K = 1,p%NBAND-1
                  I=I+1
                  IF(I.GE.1) THEN
                      J=MM-I+1
                      misc%GLF(LEG,I)=misc%GLF(LEG,I)-misc%GLK(LEG,I,J)*VP
                      misc%GLK(LEG,I,J)=0.
                  ENDIF
              ENDDO
              misc%GLK(LEG,MM,1)=1.
              misc%GLF(LEG,MM)=VP
              I=MM
              DO K = 2,p%NBAND
                  I=I+1
                  IF(I.LE.p%NEQ(LEG)) THEN
                      misc%GLF(LEG,I)=misc%GLF(LEG,I)-misc%GLK(LEG,MM,K)*VP
                      misc%GLK(LEG,MM,K)=0.
                  ENDIF
              ENDDO
          ENDDO  
       ENDDO ! End LEG loop

       ! Solver
       CALL Solve(p, x, OtherState, ErrStat, ErrMsg)

       ! 1. STATIC ITERATION (CONVERGENCE)
       IF(.NOT.DYN) THEN
           SUM1=0.0
           SUM2=0.0
           SUM3=0.0
           SUMA=0.0
           SUMB=0.0
           SUMC=0.0
           DO LEG = 1,p%NumLines
               DO I = 1,p%NEQ(LEG)
                   x%GLU(LEG,I)=x%GLU(LEG,I)+misc%GLF(LEG,I)
               ENDDO

               DO I = 1,p%NumElems+1
                   DO J = 1,p%NDIM
                       SUM1=SUM1+misc%GLF(LEG,(I-1)*p%NHBD+2*J-1)**2
                       SUMA=SUMA+x%GLU(LEG,(I-1)*p%NHBD+2*J-1)**2
                       SUM2=SUM2+misc%GLF(LEG,(I-1)*p%NHBD+2*J)**2
                       SUMB=SUMB+x%GLU(LEG,(I-1)*p%NHBD+2*J)**2
                   ENDDO
                   SUM3=SUM3+misc%GLF(LEG,I*p%NHBD-1)**2
                   SUMC=SUMC+x%GLU(LEG,I*p%NHBD-1)**2          
               ENDDO
           ENDDO

           ERR1=SQRT(SUM1/SUMA)
           ERR2=SQRT(SUM2/SUMB)
           ERR3=SQRT(SUM3/SUMC)

           IF(ERR1.GT.p%Eps.OR.ERR2.GT.p%Eps.OR.ERR3.GT.p%Eps.OR.OtherState%INCR.EQ.0) THEN
               OtherState%INCR=OtherState%INCR+1
               IF(OtherState%INCR.GT.p%MaxIter) THEN
                   ErrStat = ErrID_Fatal
                   ErrMsg = 'FEAM_Solve:No Convergence'
                   RETURN
               END IF
           ELSE IF(.NOT.OtherState%BottomTouch) THEN
               call wrscr('Number of iterations without Seabed :'//trim(num2lstr(OtherState%INCR)))
               OtherState%BottomTouch = .TRUE.
               OtherState%INCR  = 0
               Iteration1 = .TRUE.
           ELSE IF(OtherState%BottomTouch) THEN    
               call wrScr( 'Number of final iterations          :'//trim(num2lstr(OtherState%INCR)))
               OtherState%INCR = 0 
               Iteration2 = .TRUE.    

               !SET THE STATIC VALUE TO LINE MESH
               DO I = 1, p%NumLines
                   DO NN = 1,p%NumElems+1
                       DO J = 1,p%NDIM
                           MM = p%IDOF(J,1)+(NN-1)*p%NHBD
                           misc%Line_Coordinate(I,NN,J) = x%GLU(I,MM)
                           misc%Line_Tangent(I,NN,J) = x%GLU(I,MM+1)
                       END DO
                   END DO
               END DO            
           END IF
       END IF

       ! Run DO-LOOP Once for Dynamic Analysis
       IF(DYN) THEN
           Iteration1 = .TRUE.
           Iteration2 = .TRUE.
       END IF         

   ENDDO ! DO WHILE LOOP 

    ! 2. DYNAMIC INTEGRATION        
    IF(DYN) THEN      
        DO LEG=1,p%NumLines
            DO I=1,p%NEQ(LEG)
                IVERT=(I/8)*8+5
                OtherState%GLU0(LEG,I)=x%GLU(LEG,I)
                x%GLU(LEG,I)=x%GLU(LEG,I)+misc%GLF(LEG,I)
                TEMP=misc%GLF(LEG,I)*2.0/p%DT-x%GLDU(LEG,I)
                OtherState%GLDDU(LEG,I)=(TEMP-x%GLDU(LEG,I))/p%DT
                x%GLDU(LEG,I)=TEMP    
            END DO
        ENDDO   

        OtherState%INCR = OtherState%INCR + 1      
    ENDIF      
 
CONTAINS
    !...............................................................................................................................
    SUBROUTINE CheckError(ErrID,Msg)
    ! This subroutine sets the error message and level and cleans up if the error is >= AbortErrLev
    !...............................................................................................................................

    ! Passed arguments
    INTEGER(IntKi), INTENT(IN) :: ErrID       ! The error identifier (ErrStat)
    CHARACTER(*),   INTENT(IN) :: Msg         ! The error message (ErrMsg)

    INTEGER(IntKi)             :: ErrStat3    ! The error identifier (ErrStat)
    CHARACTER(ErrMsgLen)       :: ErrMsg3     ! The error message (ErrMsg)

    !............................................................................................................................
    ! Set error status/message;
    !............................................................................................................................

    IF ( ErrID /= ErrID_None ) THEN

    IF ( LEN_TRIM(ErrMsg) > 0 ) ErrMsg = TRIM(ErrMsg)//NewLine
    ErrMsg = TRIM(ErrMsg)//'FEAM_Solve:'//TRIM(Msg)
    ErrStat = MAX(ErrStat, ErrID)

    !.........................................................................................................................
    ! Clean up if we're going to return on error: close files, deallocate local arrays
    !.........................................................................................................................
    IF ( ErrStat >= AbortErrLev ) THEN
        !CALL FEAM_DestroyInputFile(InputFileData, ErrStat3, ErrMsg3 )
    END IF

    END IF

    END SUBROUTINE CheckError 
    !...............................................................................................................................
    SUBROUTINE Element(t, step, NumLinein,NumElemin, u, p, x, xd, z, OtherState, misc, ErrStat, ErrMsg)

    TYPE(FEAM_InputType),           INTENT(INOUT)  :: u           ! An initial guess for the input; input mesh must be defined
    TYPE(FEAM_ParameterType),       INTENT(IN   )  :: p           ! Parameters
    TYPE(FEAM_ContinuousStateType), INTENT(INOUT)  :: x           ! Initial continuous states
    TYPE(FEAM_DiscreteStateType),   INTENT(INOUT)  :: xd          ! discrete states
    TYPE(FEAM_ConstraintStateType), INTENT(INOUT)  :: z           ! constraint states
    TYPE(FEAM_OtherStateType),      INTENT(INOUT)  :: OtherState  ! other states
    TYPE(FEAM_MiscVarType),         INTENT(INOUT)  :: misc        ! misc/optimization variables

    INTEGER(IntKi),                 INTENT(  OUT)  :: ErrStat     ! Error status of the operation
    CHARACTER(*),                   INTENT(  OUT)  :: ErrMsg      ! Error message if ErrStat /= ErrID_None

    ! Passed variables:

    REAL(ReKi)                   :: A1,A2,A3,A4,A5                                     
    REAL(ReKi)                   :: TEMP,TEMP1,TEMP2,TEMP3                             
    REAL(ReKi)                   :: ELK(15,15),ELC(15,15)                              
    REAL(ReKi)                   :: RMAST(4,4),ADDMT(4,4,4,4),RDAMP(4,4),ADAMP(4,4,4,4)
    REAL(ReKi)                   :: W,HGHT,Y                                           
    REAL(ReKi)                   :: UU(4,4,3,3),POST(3),RDOT(3)
    REAL(ReKi)                   :: TIME,ABSV
    REAL(ReKi)                   :: WVV(3),WVA(3),ACN(3),VLN(3)
    REAL(ReKi)                   :: FTEMP(15),F0(3),F1(3),F2(3)
    REAL(ReKi)                   :: RMOD
    REAL(DbKi)                   :: t

    INTEGER(IntKi)               :: MM,NN,II                         ! Generic index
    INTEGER(IntKi)               :: I,J,K,L,M,N                      ! Generic index
    INTEGER(IntKi)               :: NumLinein, NumElemin, NumNodes    ! Generic index
    INTEGER(IntKi)               :: step



    ! Initialize variables for this routine

    ErrStat = ErrID_None
    ErrMsg  = ""

    !     SET COEFFICIENT FOR STATIC AND DYNAMIC CALCULATION 

    IF(.NOT.DYN) THEN
        A1 = 0.
        A2 = 1.
        A3 = 1.
        A4 = 0.
        A5 = 0.
    ELSE
        A1 = 4.0/(p%DT*p%DT)
        A2 = 2.
        A3 = 3.
        A4 =-1.
        A5 = 4.0/p%DT
        IF(OtherState%INCR.EQ.0) THEN
            A3=2.
            A4=0.
        END IF
    END IF

    !     SLOPE TIMES ELEML FOR THE FORMATION OF K,M MATRIX

    DO I=1,p%NDIM
        DO J=2,4,2
            misc%U(  I,J)=misc%U(  I,J)*p%Elength(NumLineIn)
            misc%U0( I,J)=misc%U0( I,J)*p%Elength(NumLineIn)
            misc%DU( I,J)=misc%DU( I,J)*p%Elength(NumLineIn)
            misc%DDU(I,J)=misc%DDU(I,J)*p%Elength(NumLineIn)
        END DO
    END DO

    !     CALCULATE U(J,M)*U(K,N) ORDER ARRAY FOR USE OF DOT PRODUCT

    DO J=1,p%NDIM
        DO M=1,4
            DO K=1,p%NDIM
                DO N=1,4
                    UU(M,N,J,K) = misc%U(J,M)*misc%U(K,N)
                END DO
            END DO
        END DO
    END DO

    !     CALCULATE MASS MATRIX,K1+K2 AND TANGENTIAL STIFFNESS KT0  

    DO I=1,p%NDIM
        DO J=1,4
            MM=p%IDOF(I,J)
            DO K=1,p%NDIM
                DO L=1,4
                    NN=p%IDOF(K,L)
                    IF(NN.GE.MM) THEN
                        TEMP1=0.
                        DO M=1,4
                            DO N=1,4
                                TEMP1=TEMP1+p%RADDM(M,N,J,L)*UU(M,N,I,K)
                            END DO
                        END DO
                        IF(I.EQ.K) THEN
                            TEMP1=p%RMASS(J,L)*p%Elength(NumLineIn)-TEMP1/p%Elength(NumLineIn)
                            misc%EMASS(MM,NN)=p%LDMassDen(NumLineIn)*TEMP1+p%LMassDen(NumLineIn)*p%RMASS(J,L)*p%Elength(NumLineIn)
                            TEMP2=0.
                            DO M=1,3
                                TEMP2 = TEMP2+p%STEN(M,J,L)*z%TSN(M)
                            END DO
                            ELK(MM,NN) = TEMP2/p%Elength(NumLineIn)
                            misc%ESTIF(MM,NN) = ELK(MM,NN)+A1*(A3* misc%EMASS(MM,NN)+A4*OtherState%EMAS0(MM,NN))*0.5
                        ELSE
                            misc%EMASS(MM,NN)=-p%LDMassDen(NumLineIn)*TEMP1/p%Elength(NumLineIn)
                            misc%ESTIF(MM,NN)= A1*(A3*misc%EMASS(MM,NN)+A4*OtherState%EMAS0(MM,NN))*0.5
                        END IF
                        misc%EMASS(NN,MM) = misc%EMASS(MM,NN)
                        misc%ESTIF(NN,MM) = misc%ESTIF(MM,NN)
                        ELK(NN,MM) = ELK(MM,NN)
                    END IF
                END DO
            END DO
        END DO
    END DO

    !     CALCULATE TANGENTIAL STIFFNESS KT1

    DO I=1,p%NDIM
        DO J=1,4
            MM=p%IDOF(I,J)
            DO K=1,3
                NN=p%JDOF(K)

                TEMP1=0.
                DO M=1,4
                    TEMP1=TEMP1 + p%STEN(K,J,M)*misc%U(I,M)
                END DO

                misc%ESTIF(MM,NN) = A2*TEMP1/p%Elength(NumLineIn)
                misc%ESTIF(NN,MM) = misc%ESTIF(MM,NN)
            END DO
        END DO
    END DO

    !     CALCULATE TANGENTIAL STIFFNESS DT1

    DO I=1,3
        MM=p%JDOF(I)
        DO J=1,3
            NN=p%JDOF(J)
            IF(NN.GE.MM) THEN

            misc%ESTIF(MM,NN) =-A2*p%PMPN(I,J)*p%Elength(NumLineIn)/p%LEAStiff(NumLineIn)
            misc%ESTIF(NN,MM) = misc%ESTIF(MM,NN)

            END IF
        END DO
    END DO


    !     CALCULATE TANGENTIAL STEFFNESS AND FORCE FROM POSSIBLE 
    !     TOUCHING OF A ELASTIC FOUNDATION WITH QUADRATIC SPRING
    !     IN VERTICAL DIRECTION

    DO I=1,p%NBAND
        misc%FORCE(I)=0.
    END DO

    IF(OtherState%BottomTouch) THEN

       II=p%NDIM
       DO N=1,p%NTRAP
           W=p%Elength(NumLineIn)/(p%NTRAP-1)
           IF(N.EQ.1.OR.N.EQ.p%NTRAP) W=0.5*W
           HGHT=0.
           DO J=1,4
               HGHT=HGHT+p%SHAPT(N,J)*misc%U(II,J)
           END DO
           Y=HGHT-p%BottmElev(NumLineIn)
           IF(Y.LT.0) THEN 
               DO L=1,4
                   MM=p%IDOF(II,L)
                   misc%FORCE(MM)=misc%FORCE(MM)+p%SHAPT(N,L)*W*p%BottmStiff(NumLineIn)*Y*Y
                   DO K=1,4
                       NN=p%IDOF(II,K)
                       IF(NN.GE.MM) THEN
                           misc%ESTIF(MM,NN) = misc%ESTIF(MM,NN)-p%SHAPT(N,L)*p%SHAPT(N,K)*W*p%BottmStiff(NumLineIn)*Y*2.0
                           misc%ESTIF(NN,MM) = misc%ESTIF(MM,NN)
                       END IF
                   END DO
               END DO
           END IF
       END DO
    END IF

    !     CALCULATE FORCE FROM CONSTANT DISTRIBUTED LOAD (GRAVITY) 

    DO I=1,p%NDIM
        DO L=1,4
            NN=p%IDOF(I,L)
            misc%FORCE(NN) = misc%FORCE(NN)+p%GRAV(I)*p%AM(L)*p%Elength(NumLineIn)
        END DO
    END DO

    !     CALCULATE HYDRDYNAMIC FORCE FROM MORISON EQUATION 
    !     USE GAUSS QUADRATURE INTEGRATION
    IF(DYN) THEN
        DO N=1,p%NGAUSS
            W=p%GAUSSW(N)*p%Elength(NumLineIn)
            DO K=1,p%NDIM
                TEMP1=0.
                TEMP2=0.
                DO J=1,4
                    TEMP1=TEMP1+p%SHAP(N,J)*misc%U(K,J)
                    TEMP2=TEMP2+p%SHAP(N,J)*misc%DU(K,J)
                END DO
                POST(K)=TEMP1
                RDOT(K)=TEMP2
            END DO
    
            NumNodes = (NumLinein-1)*(p%NumElems+1) + NumElemin + 1
            DO L=1, p%NDIM
                WVV(L) = InterpWrappedStpReal ( REAL(t, ReKi), p%WaveTime(:), p%WaveVel0(:, NumNodes, L), &
                                    misc%LastIndWave, p%NStepWave + 1 )
                WVA(L) = InterpWrappedStpReal ( REAL(t, ReKi), p%WaveTime(:), p%WaveAcc0(:, NumNodes, L), &
                                    misc%LastIndWave, p%NStepWave + 1 )
            END DO
    
            DO I=1,p%NDIM
                TEMP1=0.
                TEMP2=0.
                DO J=1,p%NDIM
                    DO K=1,4
                        DO M=1,4
                            TEMP1=TEMP1+WVA(J)*UU(K,M,J,I)*p%SHAPS(N,K)*p%SHAPS(N,M)
                            TEMP2=TEMP2+(WVV(J)-RDOT(J))*p%SHAPS(N,K)*p%SHAPS(N,M)*UU(K,M,J,I)
                        END DO
                    END DO
                END DO
                ACN(I)=WVA(I)-TEMP1/(p%Elength(NumLineIn)*p%Elength(NumLineIn))
                VLN(I)=(WVV(I)-RDOT(I))-TEMP2/(p%Elength(NumLineIn)*p%Elength(NumLineIn))
            END DO
    
            ABSV = 0.
            DO I=1,p%NDIM
                ABSV=ABSV+VLN(I)*VLN(I)
            END DO
            ABSV=SQRT(ABSV)
    
            DO I=1,p%NDIM
                DO L=1,4
                    MM=p%IDOF(I,L)
                    misc%FORCE(MM)=misc%FORCE(MM)+(p%LineCI(NumLineIn)*ACN(I)+p%LineCD(NumLineIn)*ABSV*VLN(I))*p%SHAP(N,L)*W
                END DO
            END DO
        END DO
    END IF

    !     CALCULATE THE RESIDUE FORCE Fil-K*U

    DO I=1,15
        OtherState%RSDF(I)=0.
    END DO

    DO I=1,p%NDIM
        DO J=1,4
            MM=p%IDOF(I,J)
            DO K=1,p%NDIM
                DO L=1,4
                    NN=p%IDOF(K,L)
                    OtherState%RSDF(MM)=OtherState%RSDF(MM)-A2*ELK(MM,NN)*misc%U(K,L)+A5*(A3*misc%EMASS(MM,NN)+A4*OtherState%EMAS0(MM,NN))*0.5*misc%DU(K,L)
                END DO
            END DO

            OtherState%RSDF(MM)=OtherState%RSDF(MM)+A3*misc%FORCE(MM)+A4*OtherState%FORC0(MM)

        END DO
    END DO

    !     CALCULATE THE RESIDUE AXIAL EXTENSION Gm BOTH FOR 
    !     EXTENSIBLE AND INEXTENSIBLE CASE. IF THE LIN
    !     IS IN WATER, THE HYDROSTATIC PRESSURE IS EXTRACTED FROM
    !     LAMBDA. THE INEXTENSIBLE IS A APPROXIMATION

    DO N=1,3
        MM=p%JDOF(N)
        TEMP1=0.
        DO J=1,p%NDIM
            DO K=1,4
                DO L=1,4
                    TEMP1=TEMP1+p%STEN(N,K,L)*UU(K,L,J,J)
                END DO
            END DO
        END DO
        TEMP1 = 0.5*(TEMP1/p%Elength(NumLineIn) - p%PM(N)*p%Elength(NumLineIn))
        TEMP2=0.
        DO J=1,3
            TEMP2 = TEMP2+p%PMPN(N,J)*(z%TSN(J)-z%TZER(J))
        END DO
        TEMP2=TEMP2*p%Elength(NumLineIn)
        OtherState%RSDF(MM) = -A2*(TEMP1-TEMP2/p%LEAStiff(NumLineIn))
    END DO

    !     MULTIPLY ROWS AND COLUMNS ASSOCIATED WITH SLOPES BY ELEMENT LENGTH

    DO I=1,p%NDIM
        DO J=2,4,2
            M=p%IDOF(I,J)
            DO N=1,p%NBAND
                misc%ESTIF(M,N) = misc%ESTIF(M,N)*p%Elength(NumLineIn)
                misc%ESTIF(N,M) = misc%ESTIF(N,M)*p%Elength(NumLineIn)
            END DO
            OtherState%RSDF(M) = OtherState%RSDF(M)*p%Elength(NumLineIn)
        END DO
    END DO

    RETURN
    END SUBROUTINE Element
    !...............................................................................................................................
    SUBROUTINE Solve(p, x, OtherState, ErrStat, ErrMsg)
    ! This subroutine sets the error message and level and cleans up if the error is >= AbortErrLev
    !...............................................................................................................................

    TYPE(FEAM_ParameterType),       INTENT(IN   )  :: p           ! Parameters
    TYPE(FEAM_ContinuousStateType), INTENT(INOUT)  :: x           ! Initial continuous states
    TYPE(FEAM_OtherStateType),      INTENT(INOUT)  :: OtherState  ! other states

    INTEGER(IntKi),                 INTENT(  OUT)  :: ErrStat     ! Error status of the operation
    CHARACTER(*),                   INTENT(  OUT)  :: ErrMsg      ! Error message if ErrStat /= ErrID_None

    INTEGER(IntKi), PARAMETER      :: NRM=1024
    INTEGER(IntKi), PARAMETER      :: NCM=15
    REAL(ReKi)                     :: BAND(NRM,NCM),RHS(NRM),FACTOR
    INTEGER(IntKi)                 :: MEQNS,NPIV,NPIVOT,LSTSUB,NROW,NCOL,ICOL,JCOL,IJK,JKI,NEQNS                                          

    ! Initialize variables for this routine

    ErrStat = ErrID_None
    ErrMsg  = ""

    DO LEG=1,p%NumLines
        NEQNS=p%NEQ(LEG)
        DO I=1,NEQNS
            RHS(I)=misc%GLF(LEG,I)
            DO J=1,p%NBAND
                BAND(I,J)=misc%GLK(LEG,I,J)
            END DO
        END DO

        ! Linesolver    
        MEQNS=NEQNS-1                              
       DO NPIV=1,MEQNS     
            NPIVOT=NPIV+1                                                         
            LSTSUB=NPIV+p%NBAND-1                                                     
            IF(LSTSUB.GT.NEQNS) THEN                                              
                LSTSUB=NEQNS                                                      
            ENDIF                                                                 

          DO NROW=NPIVOT,LSTSUB                                              

            NCOL=NROW-NPIV+1                                                      
            FACTOR=BAND(NPIV,NCOL)/BAND(NPIV,1)                                   
            DO NCOL=NROW,LSTSUB                                                
                ICOL=NCOL-NROW+1                                                      
                JCOL=NCOL-NPIV+1                                                      
                BAND(NROW,ICOL)=BAND(NROW,ICOL)-FACTOR*BAND(NPIV,JCOL)                
            end do
            RHS(NROW)=RHS(NROW)-FACTOR*RHS(NPIV)
         end do
      end do

        ! Back substitution                                                                                                                                       
        DO IJK=2,NEQNS                                                         
            NPIV=NEQNS-IJK+2                                                          

            RHS(NPIV)=RHS(NPIV)/BAND(NPIV,1)                                          
            LSTSUB=NPIV-p%NBAND+1                                                         
            IF(LSTSUB.LT.1) THEN                                                      
                LSTSUB=1                                                              
            ENDIF                                                                     
            NPIVOT=NPIV-1                                                             
            DO JKI=LSTSUB,NPIVOT                                                   
                NROW=NPIVOT-JKI+LSTSUB                                                    
                NCOL=NPIV-NROW+1                                                          
                FACTOR=BAND(NROW,NCOL)
                RHS(NROW)=RHS(NROW)-FACTOR*RHS(NPIV)
            end do
         end do
        RHS(1)=RHS(1)/BAND(1,1)

        DO I=1,NEQNS
            misc%GLF(LEG,I)=RHS(I)
        END DO
    END DO

    END SUBROUTINE Solve 
    !...............................................................................................................................
    SUBROUTINE Couple(p, x, OtherState, misc, ErrStat, ErrMsg)
    ! This subroutine sets the error message and level and cleans up if the error is >= AbortErrLev
    !...............................................................................................................................
    TYPE(FEAM_ParameterType),       INTENT(IN   )  :: p           ! Parameters
    TYPE(FEAM_ContinuousStateType), INTENT(INOUT)  :: x           ! Initial continuous states
    TYPE(FEAM_OtherStateType),      INTENT(INOUT)  :: OtherState  ! other states
    TYPE(FEAM_MiscVarType),         INTENT(INOUT)  :: misc        ! misc/optimization variables

    INTEGER(IntKi),                 INTENT(  OUT)  :: ErrStat     ! Error status of the operation
    CHARACTER(*),                   INTENT(  OUT)  :: ErrMsg      ! Error message if ErrStat /= ErrID_None

    REAL(ReKi)                   :: FORCE(3)
    REAL(ReKi)                   :: DELTA(3,3)
    REAL(ReKi)                   :: TEMP

    ! Initialize variables for this routine

    ErrStat = ErrID_None
    ErrMsg  = "" 

    !  DELTA(I,J) IS THE KRONECKER DELTA
   CALL Eye(DELTA, ErrStat, ErrMsg)   !bjj: this is not valid Fortran 2003:  DATA DELTA/ 1.0, 0.0, 0.0,  0.0, 1.0, 0.0, 0.0, 0.0, 1.0/

    !  CALCULATE CONTRIBUTION OF LINEAR SPRING
    !  SPRING FORCE = SLIN*(X + PLOC + THETA X PLOC - R) ! Equation for N(n)

    DO I=1,p%NDIM 
        TEMP=misc%FP(I)-misc%R(I)
        FORCE(I)=misc%SLIN(I)*TEMP
    END DO

    !  RIGHT HAND SIDE FOR R AND X EQUATIONS

    DO I=1,p%NDIM
        misc%RHSR(2*I-1) = FORCE(I)
        misc%RHSR(2*I) = 0.0
    END DO            

    !  STIFFNESS COEFFICIENTS Krr,Krx,Kxx

    DO I=1,p%NDIM
        DO J=1,p%NDIM
            misc%STIFR(2*I-1,2*J-1) = misc%SLIN(I)*DELTA(I,J)
        END DO
    END DO

    !  STIFFNESS COEFFICIENT Kr'r'

    DO I=1,p%NDIM
        DO J=1,p%NDIM
            misc%STIFR(2*I,2*J)= 0.0
        END DO
    END DO           

    END SUBROUTINE Couple    
    !...............................................................................................................................
END SUBROUTINE FEAM_Solve
!----------------------------------------------------------------------------------------------------------------------------------
SUBROUTINE FEAM_End    ( u, p, x, xd, z, OtherState, y, misc, ErrStat, ErrMsg )
! This routine is called at the end of the simulation.
!..................................................................................................................................

    TYPE(FEAM_InputType),           INTENT(INOUT)  :: u           ! System inputs
    TYPE(FEAM_ParameterType),       INTENT(INOUT)  :: p           ! Parameters
    TYPE(FEAM_ContinuousStateType), INTENT(INOUT)  :: x           ! Continuous states
    TYPE(FEAM_DiscreteStateType),   INTENT(INOUT)  :: xd          ! Discrete states
    TYPE(FEAM_ConstraintStateType), INTENT(INOUT)  :: z           ! Constraint states
    TYPE(FEAM_OtherStateType),      INTENT(INOUT)  :: OtherState  ! Other states
    TYPE(FEAM_OutputType),          INTENT(INOUT)  :: y           ! System outputs
    TYPE(FEAM_MiscVarType),         INTENT(INOUT)  :: misc        ! misc/optimization variables
    INTEGER(IntKi),                 INTENT(  OUT)  :: ErrStat     ! Error status of the operation
    CHARACTER(*),                   INTENT(  OUT)  :: ErrMsg      ! Error message if ErrStat /= ErrID_None

    ! Initialize ErrStat

    ErrStat = ErrID_None
    ErrMsg  = ""


    ! Place any last minute operations or calculations here:


    ! Close files here:



    ! Destroy the input data:

    CALL FEAM_DestroyInput( u, ErrStat, ErrMsg )


    ! Destroy the parameter data:

    CALL FEAM_DestroyParam( p, ErrStat, ErrMsg )


    ! Destroy the state data:

    CALL FEAM_DestroyContState(   x,           ErrStat, ErrMsg )
    CALL FEAM_DestroyDiscState(   xd,          ErrStat, ErrMsg )
    CALL FEAM_DestroyConstrState( z,           ErrStat, ErrMsg )
    CALL FEAM_DestroyOtherState(  OtherState,  ErrStat, ErrMsg )

    CALL FEAM_DestroyMisc( misc, ErrStat, ErrMsg )

    ! Destroy the output data:

    CALL FEAM_DestroyOutput( y, ErrStat, ErrMsg )




END SUBROUTINE FEAM_End
!----------------------------------------------------------------------------------------------------------------------------------
!> This is a loose coupling routine for solving constraint states, integrating continuous states, and updating discrete and other 
!! states. Continuous, constraint, discrete, and other states are updated to values at t + Interval.
SUBROUTINE FEAM_UpdateStates  ( t, n, Inputs, InputTimes, p, x, xd, z, OtherState, misc, ErrStat, ErrMsg )
!..................................................................................................................................

    REAL(DbKi),                      INTENT(IN   ) :: t               ! Current simulation time in seconds
    INTEGER(IntKi),                  INTENT(IN   ) :: n               ! Current step of the simulation: t = n*Interval
    TYPE(FEAM_InputType),            INTENT(INOUT) :: Inputs(:)       ! Inputs at InputTimes (intent out only for record keeping in ExtrapInterp)
    REAL(DbKi),                      INTENT(IN   ) :: InputTimes(:)   ! Times in seconds associated with Inputs
    TYPE(FEAM_ParameterType),        INTENT(IN   ) :: p               ! Parameters
    TYPE(FEAM_ContinuousStateType),  INTENT(INOUT) :: x               ! Input: Continuous states at t;
                                                                      !   Output: Continuous states at t  + Interval
    TYPE(FEAM_DiscreteStateType),    INTENT(INOUT) :: xd              ! Input: Discrete states at t;
                                                                      !   Output: Discrete states at t  + Interval
    TYPE(FEAM_ConstraintStateType),  INTENT(INOUT) :: z               ! Input: Discrete states at t;
                                                                      !   Output: Constraint states at t  + Interval
    TYPE(FEAM_OtherStateType),       INTENT(INOUT) :: OtherState      ! Input: Other states at t;
                                                                      !   Output: Other states at t  + Interval
    TYPE(FEAM_MiscVarType),          INTENT(INOUT) :: misc            ! misc/optimization variables
    INTEGER(IntKi),                  INTENT(  OUT) :: ErrStat         ! Error status of the operation
    CHARACTER(*),                    INTENT(  OUT) :: ErrMsg          ! Error message if ErrStat /= ErrID_None

    ! Local variables

!    TYPE(FEAM_ContinuousStateType)                 :: dxdt            ! Continuous state derivatives at t
    TYPE(FEAM_InputType)                           :: u               ! Inputs at t
!    TYPE(FEAM_InputFile)                           :: InputFileData  ! Data stored in the module's input file

    INTEGER(IntKi)                                 :: ErrStat2        ! Error status of the operation (occurs after initial error)
    CHARACTER(ErrMsgLen)                           :: ErrMsg2         ! Error message if ErrStat2 /= ErrID_None
    INTEGER(IntKi)                                 :: I               ! Generic loop index
    CHARACTER(*), PARAMETER                        :: RoutineName = 'FEAM_UpdateStates'

    ! Initialize ErrStat

    ErrStat = ErrID_None
    ErrMsg  = ""

    ! Get the inputs, based on the array of values sent by the glue code:

    CALL FEAM_CopyInput( Inputs(1), u, MESH_NEWCOPY, ErrStat2, ErrMsg2 )          
      call setErrStat(ErrStat2,ErrMsg2,ErrStat,ErrMsg,RoutineName)
      
    !bjj: this gives us extrapolated values at t+dt; but we should use FEAM_Input_ExtrapInterp
      
    IF ( ErrStat >= AbortErrLev ) THEN
       call cleanup()
        RETURN
    END IF    


      !bjj: as far as I can tell,
      ! misc%FAST_FP is absolute position of u%PtFairleadDisplacement at t+dt, 
      ! and OtherState%FAST_FPA is the absolute position of u%PtFairleadDisplacement at t
    DO I = 1,u%PtFairleadDisplacement%NNodes
        misc%FAST_FP(I,1) = u%PtFairleadDisplacement%Position(1,I) + u%PtFairleadDisplacement%TranslationDisp(1,I)  !X
        misc%FAST_FP(I,2) = u%PtFairleadDisplacement%Position(2,I) + u%PtFairleadDisplacement%TranslationDisp(2,I)  !Y 
        misc%FAST_FP(I,3) = u%PtFairleadDisplacement%Position(3,I) + u%PtFairleadDisplacement%TranslationDisp(3,I)  !Z
    END DO   

    CALL FEAM_Solve(t, n, u, p, x, xd, z, OtherState, misc, .true., ErrStat2, ErrMsg2 ) ! @yb: This is the time integration, correct?
      call setErrStat(ErrStat2,ErrMsg2,ErrStat,ErrMsg,RoutineName)
      
    !Save current fairlead position for next step integration   
    DO I = 1,u%PtFairleadDisplacement%NNodes
        OtherState%FAST_FPA(I,1) = misc%FAST_FP(I,1)
        OtherState%FAST_FPA(I,2) = misc%FAST_FP(I,2)
        OtherState%FAST_FPA(I,3) = misc%FAST_FP(I,3)
    END DO 

    call cleanup()
    
    return
contains
      subroutine cleanup()    
         CALL FEAM_DestroyInput(u,ErrStat2,ErrMsg2)    
      end subroutine cleanup
    
END SUBROUTINE FEAM_UpdateStates
!----------------------------------------------------------------------------------------------------------------------------------
!> Routine for computing outputs, used in both loose and tight coupling.
SUBROUTINE FEAM_CalcOutput    ( t, u, p, x, xd, z, OtherState, y, misc, ErrStat, ErrMsg )
!..................................................................................................................................

    REAL(DbKi),                     INTENT(IN   )  :: t           ! Current simulation time in seconds
    TYPE(FEAM_InputType),           INTENT(IN   )  :: u           ! Inputs at t
    TYPE(FEAM_ParameterType),       INTENT(IN   )  :: p           ! Parameters
    TYPE(FEAM_ContinuousStateType), INTENT(IN   )  :: x           ! Continuous states at t
    TYPE(FEAM_DiscreteStateType),   INTENT(IN   )  :: xd          ! Discrete states at t
    TYPE(FEAM_ConstraintStateType), INTENT(IN   )  :: z           ! Constraint states at t
    TYPE(FEAM_OtherStateType),      INTENT(IN   )  :: OtherState  ! Other states at t
    TYPE(FEAM_OutputType),          INTENT(INOUT)  :: y           ! Outputs computed at t (Input only so that mesh con-
                                                                  !   nectivity information does not have to be recalculated)
    TYPE(FEAM_MiscVarType),         INTENT(INOUT)  :: misc        ! misc/optimization variables
    INTEGER(IntKi),                 INTENT(  OUT)  :: ErrStat     ! Error status of the operation
    CHARACTER(*),                   INTENT(  OUT)  :: ErrMsg      ! Error message if ErrStat /= ErrID_None

    ! Local variables
    REAL(ReKi)                                     :: AllOuts(MaxOutPts)     ! All the available output channels
    INTEGER(IntKi)                                 :: I,J                    ! Generic loop index
!    INTEGER(IntKi)                                 :: K                      ! Blade index
!    INTEGER(IntKi)                                 :: ErrStat2
!    CHARACTER(ErrMsgLen)                           :: ErrMsg2
    INTEGER(IntKi)                                 :: LEG                    ! Generic index
    REAL(ReKi)                                     :: TEMP                   ! Temporary storage

    ! Initialize ErrStat

    ErrStat = ErrID_None
    ErrMsg  = ""
     
    DO LEG=1,p%NumLines
    
        ! Fairlead & Anchor tensions Output
        misc%FAIR_T(LEG) = x%GLU(LEG,(p%NumElems+1)*p%NHBD-1)
        misc%ANCH_T(LEG) = x%GLU(LEG,p%NHBD-1)
    
        ! Fairlead & Anchor angle output
        DO I=1, p%NDIM
            misc%FAIR_ANG(LEG,I) = x%GLU(LEG,p%NumElems*p%NHBD+I*2)
            misc%ANCH_ANG(LEG,I) = x%GLU(LEG,I*2)
        ENDDO
    
        ! Line element displacement & tangent vector output
        DO I=1, p%NumElems+1
            DO J=1, p%NDIM
                misc%Line_Coordinate(LEG,I,J) = x%GLU(LEG,(I-1)*p%NHBD+J*2-1)
                misc%Line_Tangent(   LEG,I,J) = x%GLU(LEG,(I-1)*p%NHBD+J*2  )
            ENDDO
        ENDDO
    
         ! OtherState%FAST_RP is set in UpdateStates (FEAM_Solve); seems like we should just recalculate it here:        
       TEMP = SQRT(OtherState%FAST_RP(LEG,1)**2+OtherState%FAST_RP(LEG,2)**2+OtherState%FAST_RP(LEG,3)**2)
       DO J=1,p%NDIM
           misc%F_Lines(LEG,J) = misc%FAIR_T(LEG)*(-OtherState%FAST_RP(LEG,J)/TEMP)
       END DO
    
    ENDDO 
    
    DO i = 1,y%PtFairleadLoad%NNodes
        y%PtFairleadLoad%Force(1,i) = misc%F_Lines(i,1) 
        y%PtFairleadLoad%Force(2,i) = misc%F_Lines(i,2)
        y%PtFairleadLoad%Force(3,i) = misc%F_Lines(i,3)
    END DO
   
    !Below Outputs are available
    !misc%FAIR_T                         : Fairlead tension [N]         
    !misc%FAIR_ANG(Line #, Direction)    : Fairlead angle   [rad]
    !misc%ANCH_T                         : Anchor tension   [N]
    !misc%ANCH_ANG(Line #, Direction)    : Anchor angle     [rad]
    
    DO i=1,MIN(10,p%NumLines)  ! 10 tension outputs
       AllOuts(FairT(i))   = misc%FAIR_T(i)/1000
       AllOuts(AnchT(i))   = misc%ANCH_T(i)/1000
    END DO
    
    !...............................................................................................................................   
    ! Place the selected output channels into the WriteOutput(:) array with the proper sign:
    !...............................................................................................................................   

    DO I = 1,p%NumOuts  ! Loop through all selected output channels
       y%WriteOutput(I) = p%OutParam(I)%SignM * AllOuts( p%OutParam(I)%Indx ) 
    ENDDO             ! I - All selected output channels

    RETURN
CONTAINS
    !...............................................................................................................................
    SUBROUTINE CheckError(ErrID,Msg)  
    ! This subroutine sets the error message and level and cleans up if the error is >= AbortErrLev
    !...............................................................................................................................

    ! Passed arguments
    INTEGER(IntKi), INTENT(IN) :: ErrID       ! The error identifier (ErrStat)
    CHARACTER(*),   INTENT(IN) :: Msg         ! The error message (ErrMsg)


    !............................................................................................................................
    ! Set error status/message;
    !............................................................................................................................

    IF ( ErrID /= ErrID_None ) THEN

    IF ( LEN_TRIM(ErrMsg) > 0 ) ErrMsg = TRIM(ErrMsg)//NewLine
    ErrMsg = TRIM(ErrMsg)//TRIM(Msg)
    ErrStat = MAX(ErrStat, ErrID)

    !.........................................................................................................................
    ! Clean up if we're going to return on error: close files, deallocate local arrays
    !.........................................................................................................................
    IF ( ErrStat >= AbortErrLev ) THEN
    END IF

    END IF


    END SUBROUTINE CheckError        
    !...............................................................................................................................
END SUBROUTINE FEAM_CalcOutput
!----------------------------------------------------------------------------------------------------------------------------------
SUBROUTINE FEAM_ReadInput( InputFileName, InputFileData, OutFileRoot, Default_DT, Default_grav, Default_WaterDens, ErrStat, ErrMsg )
! This subroutine reads the input file and stores all the data in the SrvD_InputFile structure.
! It does not perform data validation.
!..................................................................................................................................

    ! Passed variables

    CHARACTER(*), INTENT(IN)               :: InputFileName       ! Name of the input file
    CHARACTER(*), INTENT(IN)               :: OutFileRoot         ! The rootname of all the output files written by this routine.

    !BJJ MODIFIED HERE ONLY FOR TESTING:
    TYPE(FEAM_InputFile),   INTENT(OUT)    :: InputFileData       ! Data stored in the module's input file
    REAL(DbKi),             INTENT(IN)     :: Default_DT          ! The default DT (from glue code)    
    REAL(ReKi),             INTENT(IN)     :: Default_grav        ! The default gravity (from glue code)    
    REAL(ReKi),             INTENT(IN)     :: Default_WaterDens   ! The default water density (from glue code)    
    
    INTEGER(IntKi),       INTENT(OUT)      :: ErrStat             ! The error status code
    CHARACTER(*),         INTENT(OUT)      :: ErrMsg              ! The error message, if an error occurred

    ! local variables
 
    INTEGER(IntKi)                         :: UnEcho              ! Unit number for the echo file
    INTEGER(IntKi)                         :: ErrStat2            ! The error status code
    CHARACTER(ErrMsgLen)                   :: ErrMsg2             ! The error message, if an error occurred

    ! initialize values: 

    ErrStat = ErrID_None
    ErrMsg  = ''
    
   InputFileData%DT      = Default_DT
   InputFileData%gravity = Default_grav
   InputFileData%WtrDens = Default_WaterDens

    ! get the primary/platform input-file data

    CALL ReadPrimaryFile( InputFileName, InputFileData, OutFileRoot, UnEcho, ErrStat2, ErrMsg2 )
    CALL CheckError(ErrStat2,ErrMsg2)
    IF ( ErrStat >= AbortErrLev ) RETURN


    ! we may need to read additional files here (e.g., Bladed Interface)


    ! close any echo file that was opened

    IF ( UnEcho > 0 ) CLOSE( UnEcho )        

CONTAINS
    !...............................................................................................................................
    SUBROUTINE CheckError(ErrID,Msg)
    ! This subroutine sets the error message and level and cleans up if the error is >= AbortErrLev
    !...............................................................................................................................

    ! Passed arguments
    INTEGER(IntKi), INTENT(IN) :: ErrID       ! The error identifier (ErrStat)
    CHARACTER(*),   INTENT(IN) :: Msg         ! The error message (ErrMsg)


    !............................................................................................................................
    ! Set error status/message;
    !............................................................................................................................

    IF ( ErrID /= ErrID_None ) THEN

    IF ( LEN_TRIM(ErrMsg) > 0 ) ErrMsg = TRIM(ErrMsg)//NewLine
    ErrMsg = TRIM(ErrMsg)//TRIM(Msg)
    ErrStat = MAX(ErrStat, ErrID)

    !.........................................................................................................................
    ! Clean up if we're going to return on error: close files, deallocate local arrays
    !.........................................................................................................................
    IF ( ErrStat >= AbortErrLev ) THEN
        IF ( UnEcho > 0 ) CLOSE( UnEcho )
    END IF

    END IF


    END SUBROUTINE CheckError     

END SUBROUTINE FEAM_ReadInput
!----------------------------------------------------------------------------------------------------------------------------------
SUBROUTINE ReadPrimaryFile( InputFile, InputFileData, OutFileRoot, UnEc, ErrStat, ErrMsg )
! This routine reads in the primary FEAMooring input file and places the values it reads in the InputFileData structure.
!   It opens and prints to an echo file if requested.
!..................................................................................................................................


    IMPLICIT                        NONE

    ! Passed variables
    INTEGER(IntKi),     INTENT(OUT)     :: UnEc                                ! I/O unit for echo file. If > 0, file is open for writing.
    INTEGER(IntKi),     INTENT(OUT)     :: ErrStat                             ! Error status

    CHARACTER(*),       INTENT(IN)      :: InputFile                           ! Name of the file containing the primary input data
    CHARACTER(*),       INTENT(OUT)     :: ErrMsg                              ! Error message
    CHARACTER(*),       INTENT(IN)      :: OutFileRoot                         ! The rootname of the echo file, possibly opened in this routine

    TYPE(FEAM_InputFile), INTENT(INOUT) :: InputFileData                       ! All the data in the FEAMooring input file

    ! Local variables:
    INTEGER(IntKi)                :: I                                         ! loop counter
    INTEGER(IntKi)                :: J                                         ! loop counter
!    INTEGER(IntKi)                :: NumOuts                                   ! Number of output channel names read from the file 
    INTEGER(IntKi)                :: UnIn                                      ! Unit number for reading file

    INTEGER(IntKi)                :: ErrStat2 , IOS                            ! Temporary Error status
    LOGICAL                       :: Echo                                      ! Determines if an echo file should be written
    CHARACTER(ErrMsgLen)          :: ErrMsg2                                   ! Temporary Error message
    CHARACTER(1024)               :: PriPath                                   ! Path name of the primary file
    CHARACTER(1024)               :: FTitle                                    ! "File Title": the 2nd line of the input file, which contains a description of its contents
    CHARACTER(200)                :: Line                                      ! Temporary storage of a line from the input file (to compare with "default")


    ! Initialize some variables:
    ErrStat = ErrID_None
    ErrMsg  = ""

    UnEc = -1
    Echo = .FALSE.  
    CALL GetPath( InputFile, PriPath )     ! Input files will be relative to the path where the primary input file is located.

    CALL AllocAry( InputFileData%OutList, MaxOutPts, "FEAMooring Input File's Outlist", ErrStat2, ErrMsg2 )
    CALL CheckError( ErrStat2, ErrMsg2 )
    IF ( ErrStat >= AbortErrLev ) RETURN   


    ! Get an available unit number for the file.

    CALL GetNewUnit( UnIn, ErrStat2, ErrMsg2 )
    CALL CheckError( ErrStat2, ErrMsg2 )
    IF ( ErrStat >= AbortErrLev ) RETURN


    ! Open the Primary input file.

    CALL OpenFInpFile ( UnIn, InputFile, ErrStat2, ErrMsg2 )
    CALL CheckError( ErrStat2, ErrMsg2 )
    IF ( ErrStat >= AbortErrLev ) RETURN                

    ! Read the lines up/including to the "Echo" simulation control variable
    ! If echo is FALSE, don't write these lines to the echo file. 
    ! If Echo is TRUE, rewind and write on the second try.

    I = 1 !set the number of times we've read the file
    DO 
        !-------------------------- HEADER ---------------------------------------------

        CALL ReadCom( UnIn, InputFile, 'File header: Module Version (line 1)', ErrStat2, ErrMsg2, UnEc )
        CALL CheckError( ErrStat2, ErrMsg2 )
        IF ( ErrStat >= AbortErrLev ) RETURN

        CALL ReadStr( UnIn, InputFile, FTitle, 'FTitle', 'File Header: File Description (line 2)', ErrStat2, ErrMsg2, UnEc )
        CALL CheckError( ErrStat2, ErrMsg2 )
        IF ( ErrStat >= AbortErrLev ) RETURN


        !---------------------- SIMULATION CONTROL --------------------------------------

        CALL ReadCom( UnIn, InputFile, 'Section Header: Simulation Control', ErrStat2, ErrMsg2, UnEc )
        CALL CheckError( ErrStat2, ErrMsg2 )
        IF ( ErrStat >= AbortErrLev ) RETURN

        ! Echo - Echo input to "<RootName>.ech".

        CALL ReadVar( UnIn, InputFile, Echo, 'Echo',   'Echo switch', ErrStat2, ErrMsg2, UnEc )
        CALL CheckError( ErrStat2, ErrMsg2 )
        IF ( ErrStat >= AbortErrLev ) RETURN


        IF (.NOT. Echo .OR. I > 1) EXIT !exit this loop

        ! Otherwise, open the echo file, then rewind the input file and echo everything we've read

        I = I + 1         ! make sure we do this only once (increment counter that says how many times we've read this file)

        CALL OpenEcho ( UnEc, TRIM(OutFileRoot)//'.ech', ErrStat2, ErrMsg2, FEAM_Ver )
        CALL CheckError( ErrStat2, ErrMsg2 )
        IF ( ErrStat >= AbortErrLev ) RETURN

        IF ( UnEc > 0 )  WRITE (UnEc,'(/,A,/)')  'Data from '//TRIM(FEAM_Ver%Name)//' primary input file "'//TRIM( InputFile )//'":'

        REWIND( UnIn, IOSTAT=ErrStat2 )  
        IF (ErrStat2 /= 0_IntKi ) THEN
            CALL CheckError( ErrID_Fatal, 'Error rewinding file "'//TRIM(InputFile)//'".' )      
            IF ( ErrStat >= AbortErrLev ) RETURN
        END IF         

    END DO    

    CALL WrScr( ' Heading of the '//TRIM(FEAM_Ver%Name)//' input file: '//TRIM( FTitle ) )      


      ! DT - Requested integration time for FEAMooring (seconds):
   CALL ReadVar( UnIn, InputFile, Line, "DT", "Requested integration time for FEAMooring (seconds)", ErrStat2, ErrMsg2, UnEc)
      CALL CheckError( ErrStat2, ErrMsg2 )
      IF ( ErrStat >= AbortErrLev ) RETURN
      CALL Conv2UC( Line )
      IF ( INDEX(Line, "DEFAULT" ) /= 1 ) THEN ! If it's not "default", read this variable; otherwise use the value already stored in InputFileData%DT
         READ( Line, *, IOSTAT=IOS) InputFileData%DT
            CALL CheckIOS ( IOS, InputFile, 'DT', NumType, ErrStat2, ErrMsg2 )
            CALL CheckError(ErrStat2, ErrMsg2)
            IF ( ErrStat >= AbortErrLev ) RETURN         
      END IF
        
      
    ! NumLines - Number of Mooring lines:
    CALL ReadVar( UnIn, InputFile, InputFileData%NumLines, "NumLines", "Number of Mooring lines", ErrStat2, ErrMsg2, UnEc)
    CALL CheckError( ErrStat2, ErrMsg2 )
    IF ( ErrStat >= AbortErrLev ) RETURN

    ! NumElems - Number of elements per line:
    CALL ReadVar( UnIn, InputFile, InputFileData%NumElems, "NumElems", "Number of elements per line", ErrStat2, ErrMsg2, UnEc)
    CALL CheckError( ErrStat2, ErrMsg2 )
    IF ( ErrStat >= AbortErrLev ) RETURN
    
    ! Gravity - Gravitational constant:
    CALL ReadVar( UnIn, InputFile, Line, "Gravity", "Gravitational constant", ErrStat2, ErrMsg2, UnEc)
      CALL CheckError( ErrStat2, ErrMsg2 )
      IF ( ErrStat >= AbortErrLev ) RETURN
      CALL Conv2UC( Line )
      IF ( INDEX(Line, "DEFAULT" ) /= 1 ) THEN ! If it's not "default", read this variable; otherwise use the value already stored in InputFileData%gravity
         READ( Line, *, IOSTAT=IOS) InputFileData%Gravity
            CALL CheckIOS ( IOS, InputFile, 'Gravity', NumType, ErrStat2, ErrMsg2 )
            CALL CheckError(ErrStat2, ErrMsg2)
            IF ( ErrStat >= AbortErrLev ) RETURN                  
      END IF    
    
    
    ! WtrDens - Water density:
    CALL ReadVar( UnIn, InputFile, Line, "WtrDens", "Water density", ErrStat2, ErrMsg2, UnEc)
    CALL CheckError( ErrStat2, ErrMsg2 )
    IF ( ErrStat >= AbortErrLev ) RETURN 
      CALL Conv2UC( Line )
      IF ( INDEX(Line, "DEFAULT" ) /= 1 ) THEN ! If it's not "default", read this variable; otherwise use the value already stored in InputFileData%WtrDens
         READ( Line, *, IOSTAT=IOS) InputFileData%WtrDens
            CALL CheckIOS ( IOS, InputFile, 'WtrDens', NumType, ErrStat2, ErrMsg2 )
            CALL CheckError(ErrStat2, ErrMsg2)
            IF ( ErrStat >= AbortErrLev ) RETURN                  
      END IF  
      
    ! MaxIter - Maximum number of static iterations:
    CALL ReadVar( UnIn, InputFile, InputFileData%MaxIter, "MaxIter", "Maximum number of static iterations", ErrStat2, ErrMsg2, UnEc)
    CALL CheckError( ErrStat2, ErrMsg2 )
    IF ( ErrStat >= AbortErrLev ) RETURN 
    
    ! Eps - Static iteration tolerance:
    CALL ReadVar( UnIn, InputFile, InputFileData%Eps, "Eps", "Static iteration tolerance", ErrStat2, ErrMsg2, UnEc)
    CALL CheckError( ErrStat2, ErrMsg2 )
    IF ( ErrStat >= AbortErrLev ) RETURN 

    ! Allocate arrays for input, based on maximum allowed number of blades and outputs
    CALL AllocAry( InputFileData%LEAStiff, InputFileData%NumLines, 'Axial Stiffness array', ErrStat2, ErrMsg2 )
    CALL CheckError( ErrStat2, ErrMsg2 )
    IF ( ErrStat >= AbortErrLev ) RETURN
    
    CALL AllocAry( InputFileData%LineCI, InputFileData%NumLines, 'LineCI array', ErrStat2, ErrMsg2 )
    CALL CheckError( ErrStat2, ErrMsg2 )
    IF ( ErrStat >= AbortErrLev ) RETURN
    
    CALL AllocAry( InputFileData%LineCD, InputFileData%NumLines, 'LineCD array', ErrStat2, ErrMsg2 )
    CALL CheckError( ErrStat2, ErrMsg2 )
    IF ( ErrStat >= AbortErrLev ) RETURN

    CALL AllocAry( InputFileData%LMassDen, InputFileData%NumLines, 'Mass Density array', ErrStat2, ErrMsg2 )
    CALL CheckError( ErrStat2, ErrMsg2 )
    IF ( ErrStat >= AbortErrLev ) RETURN      

    CALL AllocAry( InputFileData%LDMassDen, InputFileData%NumLines, 'Displaced Mass Density array', ErrStat2, ErrMsg2 )
    CALL CheckError( ErrStat2, ErrMsg2 )
    IF ( ErrStat >= AbortErrLev ) RETURN

    CALL AllocAry( InputFileData%BottmStiff, InputFileData%NumLines, 'Bottom Stiffness array', ErrStat2, ErrMsg2 )
    CALL CheckError( ErrStat2, ErrMsg2 )
    IF ( ErrStat >= AbortErrLev ) RETURN

    CALL AllocAry( InputFileData%LUnstrLen, InputFileData%NumLines, 'Unstretched Length', ErrStat2, ErrMsg2 )
    CALL CheckError( ErrStat2, ErrMsg2 )
    IF ( ErrStat >= AbortErrLev ) RETURN

    CALL AllocAry( InputFileData%LRadAnch, InputFileData%NumLines, 'Anchor Radius array', ErrStat2, ErrMsg2 )
    CALL CheckError( ErrStat2, ErrMsg2 )
    IF ( ErrStat >= AbortErrLev ) RETURN

    CALL AllocAry( InputFileData%LAngAnch, InputFileData%NumLines, 'Anchor Angle array', ErrStat2, ErrMsg2 )
    CALL CheckError( ErrStat2, ErrMsg2 )
    IF ( ErrStat >= AbortErrLev ) RETURN

    CALL AllocAry( InputFileData%LDpthAnch, InputFileData%NumLines, 'Anchor Depth array', ErrStat2, ErrMsg2 )
    CALL CheckError( ErrStat2, ErrMsg2 )
    IF ( ErrStat >= AbortErrLev ) RETURN

    CALL AllocAry( InputFileData%LRadFair, InputFileData%NumLines, 'Fairlead Radius array', ErrStat2, ErrMsg2 )
    CALL CheckError( ErrStat2, ErrMsg2 )
    IF ( ErrStat >= AbortErrLev ) RETURN

    CALL AllocAry( InputFileData%LAngFair, InputFileData%NumLines, 'Fairlead Angle array', ErrStat2, ErrMsg2 )
    CALL CheckError( ErrStat2, ErrMsg2 )
    IF ( ErrStat >= AbortErrLev ) RETURN

    CALL AllocAry( InputFileData%LDrftFair, InputFileData%NumLines, 'Fairlead Draft array', ErrStat2, ErrMsg2 )
    CALL CheckError( ErrStat2, ErrMsg2 )
    IF ( ErrStat >= AbortErrLev ) RETURN

    CALL AllocAry( InputFileData%Tension, InputFileData%NumLines, 'Line Top Tension array', ErrStat2, ErrMsg2 )
    CALL CheckError( ErrStat2, ErrMsg2 )
    IF ( ErrStat >= AbortErrLev ) RETURN   

    CALL AllocAry( InputFileData%GSL, InputFileData%NumLines, 2, 3, 'Linear Stiffness array', ErrStat2, ErrMsg2 )
    CALL CheckError( ErrStat2, ErrMsg2 )
    IF ( ErrStat >= AbortErrLev ) RETURN  

    DO J=1,InputFileData%NumLines

    !---------------------- MOORING PROPERTIES -------------------------------------------
    CALL ReadCom( UnIn, InputFile, 'Section Header: Material Properties', ErrStat2, ErrMsg2, UnEc )
    CALL CheckError( ErrStat2, ErrMsg2 )
    IF ( ErrStat >= AbortErrLev ) RETURN

    ! LEAStiff - Axial Stiffness:
    CALL ReadVar( UnIn, InputFile, InputFileData%LEAStiff(J), "LEAStiff", "Axial Stiffness", ErrStat2, ErrMsg2, UnEc)
    CALL CheckError( ErrStat2, ErrMsg2 )
    IF ( ErrStat >= AbortErrLev ) RETURN

    ! LMassDen - Mass Density:
    CALL ReadVar( UnIn, InputFile, InputFileData%LMassDen(J), "LMassDen", "Mass Density", ErrStat2, ErrMsg2, UnEc)
    CALL CheckError( ErrStat2, ErrMsg2 )
    IF ( ErrStat >= AbortErrLev ) RETURN

    ! LDMassDen - Displaced Mass Density:
    CALL ReadVar( UnIn, InputFile, InputFileData%LDMassDen(J), "LDMassDen", "Displaced Mass Density", ErrStat2, ErrMsg2, UnEc)
    CALL CheckError( ErrStat2, ErrMsg2 )
    IF ( ErrStat >= AbortErrLev ) RETURN

    ! LineCI - Inertia Coefficient:
    CALL ReadVar( UnIn, InputFile, InputFileData%LineCI(J), "LineCI", "Inertia Coefficient", ErrStat2, ErrMsg2, UnEc)
    CALL CheckError( ErrStat2, ErrMsg2 )
    IF ( ErrStat >= AbortErrLev ) RETURN

    ! LineCD - Drag Coefficient:
    CALL ReadVar( UnIn, InputFile, InputFileData%LineCD(J), "LineCD", "Drag Coefficient", ErrStat2, ErrMsg2, UnEc)
    CALL CheckError( ErrStat2, ErrMsg2 )
    IF ( ErrStat >= AbortErrLev ) RETURN

    ! LUnstrLen - Total length:
    CALL ReadVar( UnIn, InputFile, InputFileData%LUnstrLen(J), "LUnstrLen", "Total length", ErrStat2, ErrMsg2, UnEc)
    CALL CheckError( ErrStat2, ErrMsg2 )
    IF ( ErrStat >= AbortErrLev ) RETURN 

    ! BottmStiff - Bottom Stiffness:
    CALL ReadVar( UnIn, InputFile, InputFileData%BottmStiff(J), "BottmStiff", "Bottom Stiffness", ErrStat2, ErrMsg2, UnEc)
    CALL CheckError( ErrStat2, ErrMsg2 )
    IF ( ErrStat >= AbortErrLev ) RETURN

    ! LRadAnch - Anchor Radius:
    CALL ReadVar( UnIn, InputFile, InputFileData%LRadAnch(J), "LRadAnch", "Anchor Radius", ErrStat2, ErrMsg2, UnEc)
    CALL CheckError( ErrStat2, ErrMsg2 )
    IF ( ErrStat >= AbortErrLev ) RETURN

    ! LAngAnch - Anchor Angle:
    CALL ReadVar( UnIn, InputFile, InputFileData%LAngAnch(J), "LRadAnch", "Anchor Angle", ErrStat2, ErrMsg2, UnEc)
    CALL CheckError( ErrStat2, ErrMsg2 )
    IF ( ErrStat >= AbortErrLev ) RETURN

    ! LDpthAnch - Anchor Depth:
    CALL ReadVar( UnIn, InputFile, InputFileData%LDpthAnch(J), "LDpthAnch", "Anchor Depth", ErrStat2, ErrMsg2, UnEc)
    CALL CheckError( ErrStat2, ErrMsg2 )
    IF ( ErrStat >= AbortErrLev ) RETURN

    ! LRadFair - Fairlead Radius:
    CALL ReadVar( UnIn, InputFile, InputFileData%LRadFair(J), "LRadFair", "Fairlead Radius", ErrStat2, ErrMsg2, UnEc)
    CALL CheckError( ErrStat2, ErrMsg2 )
    IF ( ErrStat >= AbortErrLev ) RETURN

    ! LAngFair - Fairlead Radius:
    CALL ReadVar( UnIn, InputFile, InputFileData%LAngFair(J), "LAngFair", "Fairlead Angle", ErrStat2, ErrMsg2, UnEc)
    CALL CheckError( ErrStat2, ErrMsg2 )
    IF ( ErrStat >= AbortErrLev ) RETURN

    ! LDrftFair - Fairlead Draft:
    CALL ReadVar( UnIn, InputFile, InputFileData%LDrftFair(J), "LDrftFair", "Fairlead Draft", ErrStat2, ErrMsg2, UnEc)
    CALL CheckError( ErrStat2, ErrMsg2 )
    IF ( ErrStat >= AbortErrLev ) RETURN

    ! Tension - Line Top Tension:
    CALL ReadVar( UnIn, InputFile, InputFileData%Tension(J), "Tension", "Line Top Tension", ErrStat2, ErrMsg2, UnEc)
    CALL CheckError( ErrStat2, ErrMsg2 )
    IF ( ErrStat >= AbortErrLev ) RETURN

    ! GSL - Linear Stiffness:
    CALL ReadAryLines( UnIn, InputFile, InputFileData%GSL(J,2,:), SIZE(InputFileData%GSL,3), "GSL", "Linear Stiffness", ErrStat2, ErrMsg2, UnEc)
    CALL CheckError( ErrStat2, ErrMsg2 )
    IF ( ErrStat >= AbortErrLev ) RETURN 

    InputFileData%LAngAnch(J)                           =  InputFileData%LAngAnch(J)*D2R                            ! Convert the azimuth angle of the current
    InputFileData%LAngFair(J)                           =  InputFileData%LAngFair(J)*D2R                            ! anchor and fairlead from degrees to radians.      

    ENDDO

    !---------------------- OUTPUT --------------------------------------------------         
    CALL ReadCom( UnIn, InputFile, 'Section Header: Output', ErrStat2, ErrMsg2, UnEc )
    CALL CheckError( ErrStat2, ErrMsg2 )
    IF ( ErrStat >= AbortErrLev ) RETURN

    ! SumPrint - Print summary data to <RootName>.sum (flag):
    CALL ReadVar( UnIn, InputFile, InputFileData%SumPrint, "SumPrint", "Print summary data to <RootName>.sum (flag)", ErrStat2, ErrMsg2, UnEc)
    CALL CheckError( ErrStat2, ErrMsg2 )
    IF ( ErrStat >= AbortErrLev ) RETURN

    ! OutFile - Switch to determine where output will be placed: (1: in module output file only; 2: in glue code output file only; 3: both) (-):
    CALL ReadVar( UnIn, InputFile, InputFileData%OutFile, "OutFile", "Switch to determine where output will be placed: {1: in module output file only; 2: in glue code output file only; 3: both} (-)", ErrStat2, ErrMsg2, UnEc)
    CALL CheckError( ErrStat2, ErrMsg2 )
    IF ( ErrStat >= AbortErrLev ) RETURN

    !   ! OutFileFmt - Format for module tabular (time-marching) output: (1: text file [<RootName>.out], 2: binary file [<RootName>.outb], 3: both):
    !CALL ReadVar( UnIn, InputFile, InputFileData%OutFileFmt, "OutFileFmt", "Format for module tabular (time-marching) output: (1: text file [<RootName>.out], 2: binary file [<RootName>.outb], 3: both)", ErrStat2, ErrMsg2, UnEc)
    !   CALL CheckError( ErrStat2, ErrMsg2 )
    !   IF ( ErrStat >= AbortErrLev ) RETURN      

    ! TabDelim - Flag to cause tab-delimited text output (delimited by space otherwise) (flag):
    CALL ReadVar( UnIn, InputFile, InputFileData%TabDelim, "TabDelim", "Flag to cause tab-delimited text output (delimited by space otherwise) (flag)", ErrStat2, ErrMsg2, UnEc)
    CALL CheckError( ErrStat2, ErrMsg2 )
    IF ( ErrStat >= AbortErrLev ) RETURN

    ! OutFmt - Format used for module's text tabult output (except time); resulting field should be 10 characters (-):
    CALL ReadVar( UnIn, InputFile, InputFileData%OutFmt, "OutFmt", "Format used for module's text tabular output (except time); resulting field should be 10 characters (-)", ErrStat2, ErrMsg2, UnEc)
    CALL CheckError( ErrStat2, ErrMsg2 )
    IF ( ErrStat >= AbortErrLev ) RETURN

    ! Tstart - Time to start module's tabular output (seconds):
    CALL ReadVar( UnIn, InputFile, InputFileData%Tstart, "Tstart", "Time to start module's tabular output (seconds)", ErrStat2, ErrMsg2, UnEc)
    CALL CheckError( ErrStat2, ErrMsg2 )
    IF ( ErrStat >= AbortErrLev ) RETURN
    !
    !   ! DecFact - Decimation factor for module's tabular output (1=output every step) (-):
    !CALL ReadVar( UnIn, InputFile, InputFileData%DecFact, "DecFact", "Decimation factor for module's tabular output (1=output every step) (-)", ErrStat2, ErrMsg2, UnEc)
    !   CALL CheckError( ErrStat2, ErrMsg2 )
    !   IF ( ErrStat >= AbortErrLev ) RETURN

    !---------------------- OUTLIST  --------------------------------------------
    CALL ReadCom( UnIn, InputFile, 'Section Header: OutList', ErrStat2, ErrMsg2, UnEc )
    CALL CheckError( ErrStat2, ErrMsg2 )
    IF ( ErrStat >= AbortErrLev ) RETURN

    ! OutList - List of user-requested output channels (-):
    CALL ReadOutputList ( UnIn, InputFile, InputFileData%OutList, InputFileData%NumOuts, 'OutList', "List of user-requested output channels", ErrStat2, ErrMsg2, UnEc  )     ! Routine in NWTC Subroutine Library
    CALL CheckError( ErrStat2, ErrMsg2 )
    IF ( ErrStat >= AbortErrLev ) RETURN     

    !---------------------- END OF FILE -----------------------------------------

    CLOSE ( UnIn )
    RETURN


CONTAINS
    !...............................................................................................................................
    SUBROUTINE CheckError(ErrID,Msg)
    ! This subroutine sets the error message and level
    !...............................................................................................................................

    ! Passed arguments
    INTEGER(IntKi), INTENT(IN) :: ErrID       ! The error identifier (ErrStat)
    CHARACTER(*),   INTENT(IN) :: Msg         ! The error message (ErrMsg)


    !............................................................................................................................
    ! Set error status/message;
    !............................................................................................................................

    IF ( ErrID /= ErrID_None ) THEN

    IF ( LEN_TRIM(ErrMsg) > 0 ) ErrMsg = TRIM(ErrMsg)//NewLine
    ErrMsg = TRIM(ErrMsg)//'ReadPrimaryFile:'//TRIM(Msg)
    ErrStat = MAX(ErrStat, ErrID)

    !.........................................................................................................................
    ! Clean up if we're going to return on error: close file, deallocate local arrays
    !.........................................................................................................................
    IF ( ErrStat >= AbortErrLev ) THEN
        CLOSE( UnIn )
        !            IF ( UnEc > 0 ) CLOSE ( UnEc )
    END IF

    END IF


    END SUBROUTINE CheckError
    !...............................................................................................................................
END SUBROUTINE ReadPrimaryFile      
!----------------------------------------------------------------------------------------------------------------------------------
SUBROUTINE FEAM_SetParameters( InputFileData, p, ErrStat, ErrMsg )
! This subroutine sets the parameters, based on the data stored in InputFileData
!..................................................................................................................................

    TYPE(FEAM_InputFile),     INTENT(IN)       :: InputFileData  ! Data stored in the module's input file
    TYPE(FEAM_ParameterType), INTENT(INOUT)    :: p              ! The module's parameter data
    INTEGER(IntKi),           INTENT(OUT)      :: ErrStat        ! The error status code
    CHARACTER(*),             INTENT(OUT)      :: ErrMsg         ! The error message, if an error occurred

    ! Local variables  
    INTEGER(IntKi)                             :: ErrStat2       ! Temporary error ID   
    CHARACTER(ErrMsgLen)                       :: ErrMsg2        ! Temporary message describing error

    ! Initialize variables

    ErrStat = ErrID_None
    ErrMsg  = ''

    p%DT = InputFileData%DT

    ! Set parameters from primary input file
    CALL SetPrimaryParameters( p, InputFileData, ErrStat2, ErrMsg2  )
    CALL CheckError( ErrStat2, ErrMsg2 )
    IF ( ErrStat >= AbortErrLev ) RETURN  

       !.............................................
       ! Parameters for file output
       !.............................................
    p%NumOuts = InputFileData%NumOuts
       
    CALL SetOutParam(InputFileData%OutList, p, ErrStat2, ErrMsg2 )
       CALL CheckError(ErrStat2,ErrMsg2)
       IF (ErrStat >= AbortErrLev) RETURN  
       
    IF ( InputFileData%TabDelim ) THEN
       p%Delim = TAB
    ELSE
       p%Delim = ' '
    END IF           


CONTAINS
    !...............................................................................................................................
    SUBROUTINE CheckError(ErrID,Msg)
    ! This subroutine sets the error message and level
    !...............................................................................................................................

    ! Passed arguments
    INTEGER(IntKi), INTENT(IN) :: ErrID       ! The error identifier (ErrStat)
    CHARACTER(*),   INTENT(IN) :: Msg         ! The error message (ErrMsg)


    !............................................................................................................................
    ! Set error status/message;
    !............................................................................................................................

    IF ( ErrID /= ErrID_None ) THEN

    IF ( LEN_TRIM(ErrMsg) > 0 ) ErrMsg = TRIM(ErrMsg)//NewLine
    ErrMsg = TRIM(ErrMsg)//TRIM(Msg)
    ErrStat = MAX(ErrStat, ErrID)

    !.........................................................................................................................
    ! Clean up if we're going to return on error: close files, deallocate local arrays
    !.........................................................................................................................
    IF ( ErrStat >= AbortErrLev ) THEN
    END IF

    END IF


    END SUBROUTINE CheckError

END SUBROUTINE FEAM_SetParameters  
!---------------------------------------------------------------------------------------------------------------------------------
SUBROUTINE SetPrimaryParameters( p, InputFileData, ErrStat, ErrMsg  )
! This takes the primary input file data and sets the corresponding parameters.
!..................................................................................................................................

    IMPLICIT                        NONE


    ! Passed variables

    TYPE(FEAM_ParameterType), INTENT(INOUT)  :: p                            ! Parameters of the structural dynamics module
    TYPE(FEAM_InputFile),     INTENT(IN)     :: InputFileData                ! Data stored in the module's input file
    INTEGER(IntKi),           INTENT(OUT)    :: ErrStat                      ! Error status
    CHARACTER(*),             INTENT(OUT)    :: ErrMsg                       ! Error message

    REAL(ReKi)                               :: GPS(6),SHAPSS(6,4),PS(6,3)
    REAL(ReKi)                               :: W,S,SS,SSS 
    INTEGER(IntKi)                           :: I                      ! Generic loop index
    INTEGER(IntKi)                           :: M                      ! Generic loop index
    INTEGER(IntKi)                           :: N                      ! Generic loop index
    INTEGER(IntKi)                           :: J                      ! Generic loop index
    INTEGER(IntKi)                           :: K                      ! Generic loop index
    INTEGER(IntKi)                           :: L                      ! Generic loop index
    INTEGER(IntKi)                           :: NN, LEG

    !bjj: ERROR CHECKING!!!

    ! Initialize error data
    ErrStat = ErrID_None
    ErrMsg  = ''

    !...............................................................................................................................
    ! Direct copy of scalar variables variables:
    !...............................................................................................................................
    p%NumLines = InputFileData%NumLines
    p%NumElems = InputFileData%NumElems
    p%NDIM     = 3
    p%NHBD     = 8
    
    p%Gravity  = InputFileData%Gravity       ! Gravity
    p%WtrDens  = InputFileData%WtrDens       ! Water density
    p%MaxIter  = InputFileData%MaxIter       ! Number of maximum static iteration 
    p%Eps      = InputFileData%Eps           ! Static iteration tolerance
    
        
    !...............................................................................................................................
    ! allocate arrays
    !...............................................................................................................................
        
    CALL AllocAry( p%GSL, p%NumLines, 2, 3, 'Linear Stiffness array', ErrStat, ErrMsg )
    IF ( ErrStat >= AbortErrLev ) RETURN  

    CALL AllocAry( p%GP, p%NumLines, 3, 'GP',  ErrStat, ErrMsg )
    IF (ErrStat >= AbortErrLev) RETURN

    CALL AllocAry( p%NEQ,  p%NumLines, 'NEQ',  ErrStat, ErrMsg )
    IF (ErrStat >= AbortErrLev) RETURN

    CALL AllocAry( p%Elength,  p%NumLines, 'Elength',  ErrStat, ErrMsg )
    IF (ErrStat >= AbortErrLev) RETURN

    CALL AllocAry( p%BottmElev,  p%NumLines, 'BottmElev',  ErrStat, ErrMsg )
    IF (ErrStat >= AbortErrLev) RETURN

    CALL AllocAry( p%BottmStiff,  p%NumLines, 'BottmStiff',  ErrStat, ErrMsg )
    IF (ErrStat >= AbortErrLev) RETURN

    CALL AllocAry( p%LMassDen,  p%NumLines, 'LMassDen',  ErrStat, ErrMsg )
    IF (ErrStat >= AbortErrLev) RETURN

    CALL AllocAry( p%LDMassDen,  p%NumLines, 'LDMassDen',  ErrStat, ErrMsg )
    IF (ErrStat >= AbortErrLev) RETURN

    CALL AllocAry( p%LEAStiff,  p%NumLines, 'LEAStiff',  ErrStat, ErrMsg )
    IF (ErrStat >= AbortErrLev) RETURN
    
    CALL AllocAry( p%LineCI,  p%NumLines, 'LineCI',  ErrStat, ErrMsg )
    IF (ErrStat >= AbortErrLev) RETURN
    
    CALL AllocAry( p%LineCD,  p%NumLines, 'LineCD',  ErrStat, ErrMsg )
    IF (ErrStat >= AbortErrLev) RETURN

    CALL AllocAry( p%Bvp,  p%NumLines, 3, 'Bvp',  ErrStat, ErrMsg )
    IF (ErrStat >= AbortErrLev) RETURN

    CALL AllocAry( p%GLUZR,  p%NumLines, p%NumElems+1, 2*p%NDIM, 'GLUZR',  ErrStat, ErrMsg )
    IF (ErrStat >= AbortErrLev) RETURN

    CALL AllocAry( p%GTZER,  p%NumLines, 2*p%NumElems+1, 'GTZER',  ErrStat, ErrMsg )
    IF (ErrStat >= AbortErrLev) RETURN
    
    !...............................................................................................................................
    ! Direct copy of variables:
    !...............................................................................................................................
    p%NEQ(:)   = (p%NumElems+1)*p%NHBD-1     ! number of equations
    
    DO I = 1, p%NumLines
        p%GSL        (I,:,:)  = InputFileData%GSL     (I,:,:)  !bjj: todo: Inspector flags this as uninitialized
        p%LineCI     (I)      = InputFileData%LineCI  (I)
        p%LineCD     (I)      = InputFileData%LineCD  (I)
        p%Elength    (I)      = InputFileData%LUnstrLen(I)/p%NumElems
        p%BottmElev  (I)      = -InputFileData%LDpthAnch  (I)
        p%BottmStiff (I)      = InputFileData%BottmStiff (I)
        p%LMassDen   (I)      = InputFileData%LMassDen   (I)
        p%LDMassDen  (I)      = InputFileData%LDMassDen  (I)
        p%LEAStiff   (I)      = InputFileData%LEAStiff   (I)
    ENDDO       

    
    ! set p%GTZER
    
   DO I=1,2*p%NumElems+1
      p%GTZER(:,I) = InputFileData%Tension(:)
   ENDDO   
    
   DO I = 1,p%NumLines
      NN = 2*p%NumElems+1
      DO J = 3,NN-2,2
         p%GTZER(I,J) = 0.5*(p%GTZER(I,J-1)+p%GTZER(I,J+1))
      ENDDO

      p%GTZER(I,1)  = 2.0*p%GTZER(I,2   )-p%GTZER(I,3)
      p%GTZER(I,NN) = 2.0*p%GTZER(I,NN-1)-p%GTZER(I,NN-2)       
   END DO    
    !...............
    
    ! Line Static Analysis
    DO LEG = 1, p%NumLines ! Start LEG Loop     
       DO I = 1,p%NDIM
           p%GRAV(I)=0.0
       ENDDO
       p%GRAV(p%NDIM)=-p%Gravity*(p%LMassDen(LEG)-p%LDMassDen(LEG))
    ENDDO
    
    
    ! SET UP 6 POINT GAUSSIAN QUADRATURE INTEGRATION
    ! VALUES FROM ABRAMOWITZ AND STEGUN, P. 916
    ! CHANGE INTEGRATION INTERVAL FROM (-1,1) TO (0,1)

    p%NGAUSS = 6
    GPS(4) =     0.238619186083197
    GPS(5) =     0.661209386466265
    GPS(6) =     0.932469514203152
    p%GAUSSW(4) = 0.467913934572691
    p%GAUSSW(5) = 0.360761573048139
    p%GAUSSW(6) = 0.171324492379170
    GPS(1) = -GPS(6)
    GPS(2) = -GPS(5)
    GPS(3) = -GPS(4)
    p%GAUSSW(1) = p%GAUSSW(6)
    p%GAUSSW(2) = p%GAUSSW(5)
    p%GAUSSW(3) = p%GAUSSW(4)

    DO I = 1,p%NGAUSS
        GPS(I) = (GPS(I) + 1.)/2.
        p%GAUSSW(I) = p%GAUSSW(I)/2.
    ENDDO


    ! CALCULATE SHAPE FUNCTIONS
    ! SAVE VALUES AND FIRST DERIVATIVES FOR USE IN THE PROGRAM

    DO N=1,p%NGAUSS
        S = GPS(N)
        SS = S*S
        SSS = S*SS

        p%SHAP(N,1)   = 1. - 3.*SS + 2.*SSS
        p%SHAPS(N,1)  =    - 6.*S  + 6.*SS
        SHAPSS(N,1) =    - 6.    + 12.*S

        p%SHAP(N,2)   = S  - 2.*SS + SSS
        p%SHAPS(N,2)  = 1. - 4.*S  + 3.*SS
        SHAPSS(N,2) =    - 4.    + 6.*S

        p%SHAP(N,3)   = 3.*SS - 2.*SSS
        p%SHAPS(N,3)  = 6.*S  - 6.*SS
        SHAPSS(N,3) = 6.    - 12.*S

        p%SHAP(N,4)   = -SS   + SSS
        p%SHAPS(N,4)  = -2.*S + 3.*SS
        SHAPSS(N,4) = -2.   + 6.*S

        PS(N,1) = 1. - 3.*S + 2.*SS
        PS(N,2) = 4.*S - 4.*SS
        PS(N,3) = 2.*SS - S
    ENDDO

    !  CALCULATE SHAPE FUNCTIONS FOR TRANPEZOIDE INTEGRATION
    !  SAVE VALUES FOR USE IN THE PROGRAM

    p%NTRAP = 10
    DO N=1,p%NTRAP
       S = 1.0*(N-1)/REAL((p%NTRAP-1),ReKi)
        SS = S*S
        SSS = S*SS

        p%SHAPT(N,1) = 1. - 3.*SS + 2.*SSS
        p%SHAPTS(N,1)=    - 6.*S  + 6.*SS
        p%SHAPT(N,2) = S  - 2.*SS + SSS
        p%SHAPTS(N,2)= 1. - 4.*S  + 3.*SS
        p%SHAPT(N,3) = 3.*SS - 2.*SSS
        p%SHAPTS(N,3)= 6.*S  - 6.*SS
        p%SHAPT(N,4) =   -SS + SSS
        p%SHAPTS(N,4)= -2.*S + 3.*SS
    END DO

    !  CALCULATE CONSTANT ARRAYS FOR STIFFNESS AND MASS
    !  CALCULATE THE  INTEGRAL OF P(M)*P(N) 
    !  SBEND IS FOR BENDING STIFFNESS   INT EI A''A'' DS
    !  STEN IS FOR TENSION (I.E., LAGRANGE MULTIPLIER) STIFFNESS
    !  PMPN IS FOR TANGENTIAL STIFFNESS Dmn
    !  RMASS IS MULTIPLIED BY MASS + ADDED MASS
    !  RADDM IS MULTIPLIED BY -ADDED MASS TO SUBTRACT AXIAL COMPONENT

    DO I = 1,4
        p%AM(I) = 0.
        DO L = 1,4
            p%SBEND(I,L)  = 0.
            p%STEN(1,I,L) = 0.
            p%STEN(2,I,L) = 0.
            p%STEN(3,I,L) = 0.
            p%RMASS(I,L)  = 0.

            DO M = 1,4
                DO K = 1,4
                    p%RADDM(I,L,M,K) = 0.
                ENDDO
            ENDDO
        ENDDO
    ENDDO

    DO I=1,3
        p%PM(I) = 0.
        DO L = 1,3
            p%PMPN(I,L) = 0.
            DO K = 1,4
                p%PPA(I,L,K) = 0.
            ENDDO
        ENDDO
    ENDDO

    !  INTEGRATE BY GAUSSIAN QUADRATURE

    DO N = 1,p%NGAUSS
        W = p%GAUSSW(N)

        DO I = 1,4
            p%AM(I)=p%AM(I)+p%SHAP(N,I)*W

            DO L = 1,4
                p%SBEND(I,L)  = p%SBEND(I,L) +  SHAPSS(N,I)*SHAPSS(N,L)*W
                p%STEN(1,I,L) = p%STEN(1,I,L) + PS(N,1)*p%SHAPS(N,I)*p%SHAPS(N,L)*W
                p%STEN(2,I,L) = p%STEN(2,I,L) + PS(N,2)*p%SHAPS(N,I)*p%SHAPS(N,L)*W
                p%STEN(3,I,L) = p%STEN(3,I,L) + PS(N,3)*p%SHAPS(N,I)*p%SHAPS(N,L)*W
                p%RMASS(I,L)  = p%RMASS(I,L) + p%SHAP(N,I)*p%SHAP(N,L)*W

                DO M = 1,4
                    DO K = 1,4
                        p%RADDM(I,L,M,K) = p%RADDM(I,L,M,K) + p%SHAPS(N,I)*p%SHAPS(N,L)*p%SHAP(N,M)*p%SHAP(N,K)*W
                    ENDDO
                ENDDO
            ENDDO
        ENDDO

        DO I = 1,3
            p%PM(I) = p%PM(I) + PS(N,I)*W
            DO L = 1,3
                p%PMPN(I,L) = p%PMPN(I,L) + PS(N,I)*PS(N,L)*W
                DO K = 1,4
                    p%PPA(I,L,K) = p%PPA(I,L,K) + PS(N,I)*PS(N,L)*p%SHAP(N,K)*W
                ENDDO
            ENDDO
        END DO

    ENDDO

    !  SET COUNTERS
    !
    !  NDIM IS THE NUMBER OF SPACE DIMENSIONS (2 OR 3)
    !  THERE ARE NDOF = 2*NDIM DEGREES OF FREEDOM PER NODE
    !
    !  EACH DEGREE OF FREEDOM (J,I) HAS AN ASSOCIATED DEGREE OF
    !  FREEDOM NUMBER WITHIN THE ELEMENT.  IDOF(J,I) = DOF NUMBER
    !  ASSOCIATED WITH DOF (J,I). JDOF(I) = DOF NUMBER ASSOCIATED
    !  WITH LAMBDA(I)
    !                 |1  2  7  8|
    !  FOR 2D, IDOF = |3  4  9 10|, JDOF = |5 6 11|
    !  
    !                 |1  2   9  10|
    !  FOR 3D, IDOF = |3  4  11  12|, JDOF = |7 8 15|
    !                 |5  6  13  14|
    !
    N = 0
    DO J = 1,p%NDIM
        DO I = 1,2
            N = N+1
            p%IDOF(J,I) = N
        ENDDO
    ENDDO

    N = N+2
    DO J = 1,p%NDIM
        DO I = 3,4
            N = N+1
            p%IDOF(J,I) = N
        ENDDO
    ENDDO

    p%NBAND = N+1

    IF(p%NDIM.EQ.2) THEN
        p%JDOF(1) = 5
        p%JDOF(2) = 6
        p%JDOF(3) = 11
    ELSE
        p%JDOF(1) = 7
        p%JDOF(2) = 8
        p%JDOF(3) = 15
    ENDIF    

    !   ! initialize all of the DOF parameters:
    !CALL Init_DOFparameters( InputFileData, p, ErrStat, ErrMsg ) !sets p%NDOF and p%NAug
    !   IF (ErrStat >= AbortErrLev) RETURN
    !
    !   ! Set parameters for output channels:
    !CALL SetOutParam(InputFileData%OutList, p, ErrStat, ErrMsg ) ! requires: p%NumOuts, p%NumBl, p%NBlGages, p%NTwGages; sets: p%OutParam.
    !   IF (ErrStat >= AbortErrLev) RETURN
    !
    !IF ( InputFileData%TabDelim ) THEN
    !   p%Delim = TAB
    !ELSE
    !   p%Delim = ' '
    !END IF

END SUBROUTINE SetPrimaryParameters
    
!**********************************************************************************************************************************
! NOTE: The following lines of code were generated by a Matlab script called "Write_ChckOutLst.m"
!      using the parameters listed in the "OutListParameters.xlsx" Excel file. Any changes to these 
!      lines should be modified in the Matlab script and/or Excel worksheet as necessary. 
! This code was generated by Write_ChckOutLst.m at 09-Dec-2014 14:03:37.
!----------------------------------------------------------------------------------------------------------------------------------
SUBROUTINE SetOutParam(OutList, p, ErrStat, ErrMsg )
! This routine checks to see if any requested output channel names (stored in the OutList(:)) are invalid. It returns a 
! warning if any of the channels are not available outputs from the module.
!  It assigns the settings for OutParam(:) (i.e, the index, name, and units of the output channels, WriteOutput(:)).
!  the sign is set to 0 if the channel is invalid.
! It sets assumes the value p%NumOuts has been set before this routine has been called, and it sets the values of p%OutParam here.
!..................................................................................................................................

   IMPLICIT                        NONE

      ! Passed variables

   CHARACTER(ChanLen),      INTENT(IN)     :: OutList(:)                        ! The list out user-requested outputs
   TYPE(FEAM_ParameterType),    INTENT(INOUT)  :: p                                 ! The module parameters
   INTEGER(IntKi),            INTENT(OUT)    :: ErrStat                           ! The error status code
   CHARACTER(*),              INTENT(OUT)    :: ErrMsg                            ! The error message, if an error occurred

      ! Local variables

   INTEGER                      :: I                                               ! Generic loop-counting index
!   INTEGER                      :: J                                               ! Generic loop-counting index
   INTEGER                      :: INDX                                            ! Index for valid arrays

   LOGICAL                      :: CheckOutListAgain                               ! Flag used to determine if output parameter starting with "M" is valid (or the negative of another parameter)
   LOGICAL                      :: InvalidOutput(0:MaxOutPts)                      ! This array determines if the output channel is valid for this configuration
   CHARACTER(ChanLen)         :: OutListTmp                                      ! A string to temporarily hold OutList(I)

   CHARACTER(OutStrLenM1), PARAMETER  :: ValidParamAry(20) =  (/ &                  ! This lists the names of the allowed parameters, which must be sorted alphabetically
                               "ANCHT1   ","ANCHT10  ","ANCHT2   ","ANCHT3   ","ANCHT4   ","ANCHT5   ","ANCHT6   ", &
                               "ANCHT7   ","ANCHT8   ","ANCHT9   ","FAIRT1   ","FAIRT10  ","FAIRT2   ","FAIRT3   ", &
                               "FAIRT4   ","FAIRT5   ","FAIRT6   ","FAIRT7   ","FAIRT8   ","FAIRT9   "/)
   INTEGER(IntKi), PARAMETER :: ParamIndxAry(20) =  (/ &                            ! This lists the index into AllOuts(:) of the allowed parameters ValidParamAry(:)
                                   AnchT1 ,   AnchT10 ,    AnchT2 ,    AnchT3 ,    AnchT4 ,    AnchT5 ,    AnchT6 , &
                                   AnchT7 ,    AnchT8 ,    AnchT9 ,    FairT1 ,   FairT10 ,    FairT2 ,    FairT3 , &
                                   FairT4 ,    FairT5 ,    FairT6 ,    FairT7 ,    FairT8 ,    FairT9 /)
   CHARACTER(ChanLen), PARAMETER :: ParamUnitsAry(20) =  (/ &                     ! This lists the units corresponding to the allowed parameters
                               "(kN)      ","(kN)      ","(kN)      ","(kN)      ","(kN)      ","(kN)      ","(kN)      ", &
                               "(kN)      ","(kN)      ","(kN)      ","(kN)      ","(kN)      ","(kN)      ","(kN)      ", &
                               "(kN)      ","(kN)      ","(kN)      ","(kN)      ","(kN)      ","(kN)      "/)


      ! Initialize values
   ErrStat = ErrID_None
   ErrMsg = ""
   InvalidOutput = .FALSE.


!   ..... Developer must add checking for invalid inputs here: .....

!   ................. End of validity checking .................


   !-------------------------------------------------------------------------------------------------
   ! Allocate and set index, name, and units for the output channels
   ! If a selected output channel is not available in this module, set error flag.
   !-------------------------------------------------------------------------------------------------

   ALLOCATE ( p%OutParam(0:p%NumOuts) , STAT=ErrStat )
   IF ( ErrStat /= 0_IntKi )  THEN
      ErrStat = ErrID_Fatal
      ErrMsg  = "Error allocating memory for the FEAM OutParam array."
      RETURN
   ELSE
      ErrStat = ErrID_None
   ENDIF

      ! Set index, name, and units for the time output channel:

   p%OutParam(0)%Indx  = Time
   p%OutParam(0)%Name  = "Time"    ! OutParam(0) is the time channel by default.
   p%OutParam(0)%Units = "(s)"
   p%OutParam(0)%SignM = 1


      ! Set index, name, and units for all of the output channels.
      ! If a selected output channel is not available by this module set ErrStat = ErrID_Warn.

   DO I = 1,p%NumOuts

      p%OutParam(I)%Name  = OutList(I)
      OutListTmp          = OutList(I)

      ! Reverse the sign (+/-) of the output channel if the user prefixed the
      !   channel name with a "-", "_", "m", or "M" character indicating "minus".


      CheckOutListAgain = .FALSE.

      IF      ( INDEX( "-_", OutListTmp(1:1) ) > 0 ) THEN
         p%OutParam(I)%SignM = -1                         ! ex, "-TipDxc1" causes the sign of TipDxc1 to be switched.
         OutListTmp          = OutListTmp(2:)
      ELSE IF ( INDEX( "mM", OutListTmp(1:1) ) > 0 ) THEN ! We'll assume this is a variable name for now, (if not, we will check later if OutListTmp(2:) is also a variable name)
         CheckOutListAgain   = .TRUE.
         p%OutParam(I)%SignM = 1
      ELSE
         p%OutParam(I)%SignM = 1
      END IF

      CALL Conv2UC( OutListTmp )    ! Convert OutListTmp to upper case


      Indx = IndexCharAry( OutListTmp(1:OutStrLenM1), ValidParamAry )


         ! If it started with an "M" (CheckOutListAgain) we didn't find the value in our list (Indx < 1)

      IF ( CheckOutListAgain .AND. Indx < 1 ) THEN    ! Let's assume that "M" really meant "minus" and then test again
         p%OutParam(I)%SignM = -1                     ! ex, "MTipDxc1" causes the sign of TipDxc1 to be switched.
         OutListTmp          = OutListTmp(2:)

         Indx = IndexCharAry( OutListTmp(1:OutStrLenM1), ValidParamAry )
      END IF


      IF ( Indx > 0 ) THEN ! we found the channel name
         p%OutParam(I)%Indx     = ParamIndxAry(Indx)
         IF ( InvalidOutput( ParamIndxAry(Indx) ) ) THEN  ! but, it isn't valid for these settings
            p%OutParam(I)%Units = "INVALID"
            p%OutParam(I)%SignM = 0
         ELSE
            p%OutParam(I)%Units = ParamUnitsAry(Indx) ! it's a valid output
         END IF
      ELSE ! this channel isn't valid
         p%OutParam(I)%Indx  = Time                 ! pick any valid channel (I just picked "Time" here because it's universal)
         p%OutParam(I)%Units = "INVALID"
         p%OutParam(I)%SignM = 0                    ! multiply all results by zero

         ErrStat = ErrID_Warn
         ErrMsg  = p%OutParam(I)%Name//" is not an available output channel. "//TRIM(ErrMsg)
      END IF

   END DO

   RETURN
END SUBROUTINE SetOutParam
!----------------------------------------------------------------------------------------------------------------------------------
!End of code generated by Matlab script
!**********************************************************************************************************************************
!----------------------------------------------------------------------------------------------------------------------------------
SUBROUTINE CreateLineMesh( u, p, misc, ErrStat, ErrMsg )

    TYPE(FEAM_InputType),                  INTENT(INOUT)    ::  u           ! Inputs to be allocated
    TYPE(FEAM_ParameterType),              INTENT(INOUT)    ::  p           ! Parameters
    TYPE(FEAM_MiscVarType),                INTENT(INOUT)    ::  misc        ! misc/optimization variables

    INTEGER,                                INTENT(   OUT )  ::  ErrStat              ! returns a non-zero value when an error occurs  
    CHARACTER(*),                           INTENT(   OUT )  ::  ErrMsg               ! Error message if ErrStat /= ErrID_None


    REAL(R8Ki), DIMENSION(3,3)       ::  R             ! Rotation matrix transfering local frame coordinates to the global frame. @todo: do you prefer 0 index?
    REAL(ReKi), DIMENSION(3,3)       ::  USkewed       ! Skew symmetrix matrix. @todo: do you prefer 0 index?
    REAL(ReKi), DIMENSION(3,3)       ::  USkewedSquare ! Square of skew symmetrix matrix
    REAL(ReKi), DIMENSION(3)         ::  uHat          ! unit vector
    REAL(ReKi), DIMENSION(3)         ::  refAxis       ! reference axis unit vector. Should be oriented with the vertical Z global axis (opposite of gravity)
    REAL(ReKi), DIMENSION(3)         ::  lineAxis      ! vector of a straitgh-line segment parallel to the local z axis of the line element. Equals OtherState%Line_Tangent
    REAL(ReKi), DIMENSION(3)         ::  lineCrossRef  ! cross product, lineAxis (cross) refAxis 
    REAL(ReKi)                       ::  lineDotRef    ! dot product, lineAXis (dot) refAxis
    REAL(ReKi)                       ::  theta
    REAL(ReKi)                       ::  norm
    INTEGER                    ::  I, J, K

    INTEGER                    ::  numNodes 
!    REAL                       ::  Orient(3,3)

    ! Initialize ErrStat

    ErrStat = ErrID_None         
    ErrMsg  = "" 

    
    ! initialize R to zero
    !R(1,1)=0.0 , R(1,2)=0.0 , R(1,3)=0.0  
    !R(2,1)=0.0 , R(2,2)=0.0 , R(2,3)=0.0  
    !R(3,1)=0.0 , R(3,2)=0.0 , R(3,3)=0.0  
    R = 0.0

    !! initialize USkewed to zero
    !USkewed(1,1)=0.0 , USkewed(1,2)=0.0 , USkewed(1,3)=0.0  
    !USkewed(2,1)=0.0 , USkewed(2,2)=0.0 , USkewed(2,3)=0.0  
    !USkewed(3,1)=0.0 , USkewed(3,2)=0.0 , USkewed(3,3)=0.0  
    USkewed = 0.0

    ! initialize USkewed to zero
    !USkewedSquare(1,1)=0.0 , USkewedSquare(1,2)=0.0 , USkewedSquare(1,3)=0.0  
    !USkewedSquare(2,1)=0.0 , USkewedSquare(2,2)=0.0 , USkewedSquare(2,3)=0.0  
    !USkewedSquare(3,1)=0.0 , USkewedSquare(3,2)=0.0 , USkewedSquare(3,3)=0.0  
    USkewedSquare = 0.0

    ! Initialize the reference axis. This should be fixed. This is assumed to be fixed relative to the global reference frame used in FAST:
    ! the FAST inertia frame coordinate system--X positive downwind, Y to left when viewing downwind and Z vertical, opposite gravity.
    ! We want to determine the sequence of rotations needed to convert the misc%Line_Tangent
    refAxis(1)=0.0
    refAxis(2)=0.0
    refAxis(3)=1.0

    ! initialize lineAxis
    lineAxis(1)=0.0
    lineAxis(2)=0.0
    lineAxis(3)=1.0

    ! initialize cross product variable
    lineCrossRef(1)=0.0
    lineCrossRef(2)=0.0
    lineCrossRef(3)=0.0

    ! initialize uHat
    uHat(1)=0.0
    uHat(2)=0.0
    uHat(3)=1.0

    theta      = 0.0
    norm       = 0.0

           
    ! Create the input and output meshes associated with lumped loads

    p%NumNodes = p%NumLines * (p%NumElems+1)
    CALL MeshCreate( BlankMesh      = u%HydroForceLineMesh   &
    ,IOS          = COMPONENT_INPUT        &
    ,Nnodes       = p%NumNodes             &
    ,Force        = .TRUE.                 &
    ,Moment       = .TRUE.                 &
    ,ErrStat      = ErrStat                &
    ,ErrMess      = ErrMsg                 )   

    IF (ErrStat >= AbortErrLev) RETURN

    DO I=1, p%NumLines                                  ! cycle through number of lines
        DO J=1, p%NumElems+1                            ! cycle through number of element in a line

           ! START Rotation Matrix Calculation   <-------------------------------------------------------------------------+
                                                                                                                    !      |
           ! We are basically implementing this matlab pseudocode:                                                  !      |
           !                                                                                                        !      |
           !     line_axis       = [ 0.0, 1.0, 0.0 ];                                                               !      |
           !     reference_axis  = [ 1.0, 0.0, 0.0 ];                                                               !      |
           !     theta = acos( dot(line_axis,reference_axis) );                                                     !      |
           !     if norm( cross(line_axis,reference_axis) )==0                                                      !      |
           !         u = [0,0,0];                                                                                   !      |
           !     else                                                                                               !      |
           !         u = cross(line_axis,reference_axis)/norm(cross(line_axis,reference_axis));                     !      |
           !     end                                                                                                !      |
           !     u_skewed = [ 0    -u(3)  u(2) ;                                                                    !      |
           !                  u(3)  0    -u(1) ;                                                                    !      |
           !                 -u(2)  u(1)  0   ];                                                                    !      |
           !     R = eye(3) + sin(theta)*u_skewed + (1-cos(theta))*u_skewed*u_skewed                                !      |
                                                                                                                    !      |
           DO K=1, p%NDIM                                                                                           !      |
              lineAxis(K) = misc%Line_Tangent(I,J,K)                                                          !      |
           END DO                                                                                                   !      |
                                                                                                                    !      |
           lineDotRef = lineAxis(3) ! refAxis(1)*lineAxis(1) = 0, refAxis(2)*lineAxis(2) = 0,                       !      |
           ! make sure we don't have numerical issues that make the value outside +/-1
           lineDotRef = min(  1.0_ReKi, lineDotRef )
           lineDotRef = max( -1.0_ReKi, lineDotRef )
           theta = ACOS( lineDotRef )               
           
           lineCrossRef(1) = lineAxis(2)*refAXis(3) - lineAxis(3)*refAxis(2)                                        !      |
           lineCrossRef(2) = lineAxis(3)*refAXis(1) - lineAxis(1)*refAxis(3)                                        !      |
           lineCrossRef(3) = lineAxis(1)*refAXis(2) - lineAxis(2)*refAxis(1)                                        !      |
           norm = SQRT( lineCrossRef(1)**2 + lineCrossRef(2)**2 + lineCrossRef(3)**2 )                              !      |
                                                                                                                    !      |
           IF( norm.EQ.0.0 ) THEN  ! Shouldn't test for equality among REAL's... this is bad programming on my part.!      | 
              uHat(1) = 0.0                                                                                         !      |
              uHat(2) = 0.0                                                                                         !      |
              uHat(3) = 0.0                                                                                         !      |
           ELSE                                                                                                     !      |
              uHat(1) = lineCrossRef(1)/norm                                                                        !      |
              uHat(2) = lineCrossRef(2)/norm                                                                        !      |
              uHat(3) = lineCrossRef(3)/norm                                                                        !      |
           ENDIF                                                                                                    !      |
                                                                                                                    !      |
           USkewed(1,1) =  0.0
           USkewed(1,2) = -uHat(3)
           USkewed(1,3) =  uHat(2)                              !      |
           USkewed(2,1) =  uHat(3)
           USkewed(2,2) =  0.0
           USkewed(2,3) = -uHat(1)
           USkewed(3,1) = -uHat(2)
           USkewed(3,2) =  uHat(1)
           USkewed(3,3) =  0.0                                  !      |
                                                                                                                    !      |
           ! now let's square USkewed                                                                               !      |
           ! first row                                                                                              !      |
           USkewedSquare(1,1) = USkewed(1,1)*USkewed(1,1) + USkewed(1,2)*USkewed(2,1) + USkewed(1,3)*USkewed(3,1)   !      |
           USkewedSquare(1,2) = USkewed(1,1)*USkewed(1,2) + USkewed(1,2)*USkewed(2,2) + USkewed(1,3)*USkewed(3,2)   !      |
           USkewedSquare(1,3) = USkewed(1,1)*USkewed(1,3) + USkewed(1,2)*USkewed(2,3) + USkewed(1,3)*USkewed(3,3)   !      |
                                                                                                                    !      |
           ! second row                                                                                             !      |
           USkewedSquare(2,1) = USkewed(2,1)*USkewed(1,1) + USkewed(2,2)*USkewed(2,1) + USkewed(2,3)*USkewed(3,1)   !      |
           USkewedSquare(2,2) = USkewed(2,1)*USkewed(1,2) + USkewed(2,2)*USkewed(2,2) + USkewed(2,3)*USkewed(3,2)   !      |
           USkewedSquare(2,3) = USkewed(2,1)*USkewed(1,3) + USkewed(2,2)*USkewed(2,3) + USkewed(2,3)*USkewed(3,3)   !      |
                                                                                                                    !      |
           ! third row                                                                                              !      |
           USkewedSquare(3,1) = USkewed(3,1)*USkewed(1,1) + USkewed(3,2)*USkewed(2,1) + USkewed(3,3)*USkewed(3,1)   !      |
           USkewedSquare(3,2) = USkewed(3,1)*USkewed(1,2) + USkewed(3,2)*USkewed(2,2) + USkewed(3,3)*USkewed(3,2)   !      |
           USkewedSquare(3,3) = USkewed(3,1)*USkewed(1,3) + USkewed(3,2)*USkewed(2,3) + USkewed(3,3)*USkewed(3,3)   !      |
                                                                                                                    !      |
           ! lastly, calculate this:  R = eye(3) + sin(theta)*u_skewed + (1-cos(theta))*u_skewed*u_skewed           !      |
           R(1,1) = 1.0 + SIN(theta)*USkewed(1,1) + (1-COS(theta))*USkewedSquare(1,1)                                !      |
           R(1,2) = 0.0 + SIN(theta)*USkewed(1,2) + (1-COS(theta))*USkewedSquare(1,2)                                !      |
           R(1,3) = 0.0 + SIN(theta)*USkewed(1,3) + (1-COS(theta))*USkewedSquare(1,3)                                !      |
                                                                                                                    !      |
           R(2,1) = 0.0 + SIN(theta)*USkewed(2,1) + (1-COS(theta))*USkewedSquare(2,1)                                !      |
           R(2,2) = 1.0 + SIN(theta)*USkewed(2,2) + (1-COS(theta))*USkewedSquare(2,2)                                !      |
           R(2,3) = 0.0 + SIN(theta)*USkewed(2,3) + (1-COS(theta))*USkewedSquare(2,3)                                !      |
                                                                                                                     !      |
           R(3,1) = 0.0 + SIN(theta)*USkewed(3,1) + (1-COS(theta))*USkewedSquare(3,1)                                !      |
           R(3,2) = 0.0 + SIN(theta)*USkewed(3,2) + (1-COS(theta))*USkewedSquare(3,2)                                !      |
           R(3,3) = 1.0 + SIN(theta)*USkewed(3,3) + (1-COS(theta))*USkewedSquare(3,3)                                !      |
           ! END CALCULATE Rotation Matrix +-------------------------------------------------------------------------------+
           ! @todo: the reference orienttion R needs to be commited to the Mesh. 


           NumNodes = (I-1)*(p%NumElems+1) + J
           CALL MeshPositionNode (u%HydroForceLineMesh         &
                , numNodes                      &
                , misc%Line_Coordinate(I,J,:)   &
                , ErrStat                       &
                , ErrMsg                        &
                , orient=R)
           IF (ErrStat >= AbortErrLev) RETURN
           
           ! Create the mesh element
           IF ( J < p%NumElems+1) THEN
              NumNodes = (I-1)*(p%NumElems+1) + J + 1
              CALL MeshConstructElement (u%HydroForceLineMesh               &
                   , Xelement = ELEMENT_LINE2         &  
                   , P1       = numNodes - 1          &   ! node1 number
                   , P2       = numNodes              &   ! node2 number
                   , ErrStat  = ErrStat               &
                   , ErrMess  = ErrMsg                )
              IF (ErrStat >= AbortErrLev) RETURN
           END IF
           
        END DO
    
        ! add the other extra element, connecting the first node on the mooring line:
    
    
    END DO

    CALL MeshCommit ( u%HydroForceLineMesh, ErrStat, ErrMsg )   
    IF (ErrStat >= AbortErrLev) RETURN

END SUBROUTINE CreateLineMesh
!----------------------------------------------------------------------------------------------------------------------------------
SUBROUTINE CreatePtFairleadMesh( u, p, y, ErrStat, ErrMsg )

    TYPE(FEAM_InputType),                  INTENT(INOUT)    ::  u           ! Inputs to be allocated
    TYPE(FEAM_ParameterType),              INTENT(IN   )    ::  p           ! Parameters
    TYPE(FEAM_OutputType ),                INTENT(  OUT)    ::  y           ! INTENT(  OUT) : Initial system outputs (outputs are not calculated; only the output mesh is initialized)
    INTEGER,                               INTENT(   OUT )  ::  ErrStat     ! returns a non-zero value when an error occurs  
    CHARACTER(*),                          INTENT(   OUT )  ::  ErrMsg      ! Error message if ErrStat /= ErrID_None


    INTEGER                    ::  I ! J
    INTEGER                    ::  numNodes 
    REAL(ReKi)                 ::  Pos(3)

    ! Initialize ErrStat

    ErrStat = ErrID_None         
    ErrMsg  = "" 

    numNodes = p%numLines

    ! Create the input and output meshes associated with lumped loads

    CALL MeshCreate( BlankMesh         = u%PtFairleadDisplacement  &
    ,IOS             = COMPONENT_INPUT           &
    ,Nnodes          = numNodes                  &
    ,TranslationDisp = .TRUE.                    &
    ,ErrStat         = ErrStat                   &
    ,ErrMess         = ErrMsg                    )   

    IF (ErrStat >= AbortErrLev) RETURN

    ! Create the node on the mesh
    DO I=1, numNodes

       Pos(1) = p%GP(I,1)
       Pos(2) = p%GP(I,2)
       Pos(3) = p%GP(I,3)       

       CALL MeshPositionNode ( u%PtFairleadDisplacement                 &
                               , I                                      &
                               , Pos                                    &
                               , ErrStat                                &
                               , ErrMsg                                 )
       IF (ErrStat >= AbortErrLev) RETURN

       ! Create the mesh element

       CALL MeshConstructElement ( u%PtFairleadDisplacement &
                                  , ELEMENT_POINT           &
                                  , ErrStat                 &
                                  , ErrMsg                  &
                                  , I                       )
       IF (ErrStat >= AbortErrLev) RETURN

    ENDDO

    CALL MeshCommit ( u%PtFairleadDisplacement, ErrStat, ErrMsg )   
    IF (ErrStat >= AbortErrLev) RETURN  

    ! now, copy the input PtFairleadDisplacement to output            
    ! PtFairleadLoad to complete this                                 
    CALL MeshCopy ( SrcMesh  = u%PtFairleadDisplacement , &           
                    DestMesh = y%PtFairleadLoad         , &           
                    CtrlCode = MESH_SIBLING             , &           
                    IOS      = COMPONENT_OUTPUT         , &           
                    Force    = .TRUE.                   , &           
                    ErrStat  = ErrStat                  , &           
                    ErrMess  = ErrMsg                     )           

    ! End mesh initialization    


END SUBROUTINE CreatePtFairleadMesh
!----------------------------------------------------------------------------------------------------------------------------------

END MODULE FEAMooring
!**********************************************************************************************************************************
